C ********* PICKAUTOK **********************************************************
C    derived from PICKYCOR  much modified program for e.d. spot integration.
C	remember to change version number if you update.
C       VX1.00	pre-1982	JMB/RH	original version
C	VX2.00	21.11.85	TAC
C	VX2.01	18.5.87		RH	annotations added post TAC.
C	VX2.02	20.5.87		RH	larger index range up to +/-70.
C	VX3.00	25.3.92		RH	reads more lattice params from header.
C	VX3.01	08.4.92		RH	radii in Angstroms, rather than pixels.
C	VX3.02	22.5.92		RH	resolution statistics
C	VX3.03	21.3.93		RH	default NSTEP=4
C	VX3.04	12.11.94	RH	test IOVER,IUNDER in CENTRE subroutine
C	VX3.05	13.11.94	RH	change ISIZE->NX,NY and >=2 backgrounds
C	VX3.06	25.7.95 	RH	minor debug - no path to statement
C	VX3.07	11.10.95 	RH	"is OK" added to Friedel output list.
C	VX3.08	19.3.96 	RH	introduce MAXDIM for pattern size
C	VX3.09	9.05.96 	RH	consolidate Jude's improvements
C	VX3.10	6.05.97 	RH	change PEN(2) to PEN(1)
C	VX4.00	18.08.00 	RH	convert to plot2000 plot subroutines
C          "    13.6.01         TSH     P2K_FONT needed string terminator
C
C     FLAG TO SUPRESS EXCESSIVE PRINTOUT ADDED, TAC 21-NOV-85
C     NEW VAR TO DEFINE TILT ANGLE DIRECTION ON MICROSCOPE, TAC 19-NOV-84
C    NOTE:
C ****CONVENTIONS FOR TILTDIR APPLICABLE FOR 600 TO 1000 MM CAMERA 
C ****LENGTH ON EM400 MICROSCOPE, WITH OBJECTIVE LENS CURRENT UNDERFOCUSSED.
C
C     NOW DOES A YCORRECTION TO REMOVE DENSITOMETER OD DRIFT (JMB) APPROX 1984.
C     WITH ADDITIONS TO DEAL AUTOMATICALLY W. TLTAXA,TLTANG 25.6.82.
C     READS RADIAL AVERAGE BACKGROUND CURVE WITH STANDARD DEVIATIONS
C     TESTS POINTS IN BACKGROUND RASTERS AGAINST MEAN AND STNDEV AT
C     APPROPRIATE RADIUS; REPLACES VALUE BY MEAN VALUE AT THIS RADIUS 
C     IF DEVIATION IS MORE THAN 3*STNDEV; KEEPS COUNT OF NUMBER OF
C     SUCH POINTS -- MAY 82 .JMB.
C     NOW PLOTS RESIDUAL SPOT C.OF G. VECTOR ON UNIT 8, 14.1.81
C     OFFSET R**3 DISTORTION CORRECTION AND AUTOMATIC MAXRAD CUTOFF 14.1.81
C     B3 NOW SEARCHED WITHOUT OFFSET -MAY 82- MORE STABLE ALGORITHM.
C     PICKTILT WITH RADIAL BACKGROUND CORRECTION, 19.5.80
C     PICKOFF, JOYCE-LOEBL AND NIKON VERSION 29.11.79
C     LARGER DIMENSIONS AND RASTER, GREATER OVERLAP BETWEEN STRIPS, 3.12.79
C     MODIFIED FOR TILTED PATTERNS WITH BLURRED SPOTS IN ONE DIRN, 14.12.79
C     CORRECTS FOR CURVED EWALD SPHERE ROUGHLY, 20.1.81
C     CHANGED TO ADD *** TO LAST COLUMNS OF TITLE AND SOME EXTRA NUMBERS.
C     PROGRAM TO INTEGRATE SPOTS ON A LATTICE, STARTING FROM A RECTANGULAR
C     RASTER SCAN OF THE AREA . EG . OUTPUT FROM ROYAL OBSERVATORY SCANNER
C
C    DATA CARDS ARE
C    1.  A,B,G,ANGDIS,SHRINK,TILTDIR,KV,LPRINT
C         REAL SPACE CELL DIMENSIONS AND GAMMA ANGLE, ANGLE ON FILM
C         PERPENDICULAR TO WHICH A SHRINKAGE BY FACTOR SHRINK IS TO BE
C         APPLIED BEFORE CALCULATING FINAL TLTAXIS AND TILTANGLE.
C	  TILTDIR - DIRECTION OF TILT ON MICROSCOPE, 
C		    -1 = CLOCKWISE
C		    +1 = COUNTERCLOCKWISE
C	  LPRINT  - T=PRINT ALL REFLECTION INFORMATION
C    2.  J,K 
C         LOOK-UP TABLE FOR OPTICAL DENSITY; LOOKUP(J)=K
C         THIS IS LINEARLY INTERPOLATED FOR J=1,1500.  THE FIRST CARD MUST
C         HAVE J=1 AND THE LAST J=1500.
C    3.  TITLE FOR LOOK-UP TABLE  
C    4.  NPLATE
C         PLATE NUMBER            
C    5.  TITLE
C          TITLE OF DIFFRACTION PATTERN, FOR USE ON OUTPUT FILE.
C    6.  B1, B3  
C         RADIAL DISTORTION PARAMETERS, ONLY B3 USED CURRENTLY.
C    7.  X0,Y0,
C    7a. TLTAXA,TLTANG  
C         COORDINATES OF CENTRE, POSITION OF TILTAXIS (ANGLE FROM TILTAXIS
C         TO A-AXIS IN DIRECTION A TO B POSITIVE), SIZE OF TILTANGLE,
C         TLTAXA,TLTANG OVERRIDE VALUES CALCULATED INTERNALLY, UNLESS = 0.
C         NOTE-----TLTAXA HERE IS ANGLE BETWEEN TILTAXIS AND ASTAR ON FILM.
C    8.  DX1,DY1,DX2,DY2  
C        POSITION OF (1,0) AND (0,1) RELATIVE TO CENTRE. IF ZERO, THESE
C        ARE TAKEN FROM INPUT FILE HEADER. 
C    9.  ROUT,RIN,PRPMAX  : 4-Apr-1992 now in Angstroms, converted to pixels.
C         OUTER AND INNER RADII ; MAXIMUM PERPENDICULAR RESOLUTION FROM
C         TILT AXIS FOR SPOTS TO BE INTEGRATED
C    10. NXM,NYM,NSM,NXB,NYB,NSB,NPOS (TWICE),NTYPE  
C        NXMT,NYMT, ETC.
C         PEAK AND BACKGROUND RASTERS AND POSITION OF BACKGROUND.
C         THESE NUMBERS PUT IN TWICE IF IT IS A TILTED DIFFRACTION
C         PATTERN, TO GIVE RASTER SIZE, ETC AT MOST BLURRED PART OF PATTERN.
C         NPOS=1 IS AT CENTROID OF TRIANGLE, NPOS=2 IS TWICE AS NEAR SPOTS, ETC.
C         NTYPE=0,OR 1 DETERMINES POSITION OF 6 BCKGROUNDS AS SHOWN IN
C         OUTPUT OF BCKGROUND PROGRAM
C         NTYPE=2 GIVES 4 BCKGROUNDS AT CENTRES OF LATTICE; 
C    11.  FRACT,ABSOL,XAMINE,NCYC,NCYC1  
C         REFLECTIONS WITH FRIEDEL DIFFERENCES GREATER THAN FRACT AND ABSOL
C         ARE REJECTED.  WEAK REFLECTIONS ARE REJECTED FROM CENTRE OF GRAV.
C         AND LATTICE PARAMETER CALCULATIONS BY A CRITERION USING XAMINE.
C      NCYC=N FOR FIXED NUMBER OF CYCLES=N
C      NCYC=0 MEANS GO STRAIGHT TO DATA ASSESSMENT USING INPUT PARAMETERS
C      NCYC1=1 FOR PRINT ALL CENTERING RASTERS ON FIRST CYCLE
C    12.  NELIM
C         IF +VE, NUMBER OF REFLECTIONS TO BE EXCLUDED FROM LATTICE 
C         REFINEMENT; IF -VE, -NUMBER OF REFLECTIONS WHOSE RASTERS
C         ARE TO BE DISPLAYED
C    13.  NELH(I),NELK(I),I=1,8
C         UP TO 8 INDICES OF REFLECTIONS TO BE EXCLUDED FROM REFINE; OR
C         UP TO 8 INDICES OF REFLECTIONS TO BE DISPLAYED;
C         OMIT THIS CARD IF NELIM=0
C
C INPUT AND OUTPUT**************************************************************
C
C  UNIT 1 : INPUT DENSITOMETERED FILM ARRAY, any size 
C  UNIT 2 : OUTPUT OF INTEGRATED BACKGROUND SUBTRACTED INTENSITIES.
C  UNIT 3 : INPUT RADIAL DENSITY CURVE FOR BACKGROUND CORRECTION.
C  UNIT 4 : INPUT Y-AXIS DENSITOMETER DRIFT CURVE FOR BACKGROUND CORRECTION.
C  UNIT 5 : INPUT DATASTREAM.
C  UNIT 6 : OUTPUT DATASTREAM.
C  UNIT 8 : PLOTTER OUTPUT OF SPOT POSNS. AND C.OF G. DEVIATIONS.
C  UNIT 9 : INPUT header of image file with autoindexed lattice parameters
C
C
C*******************************************************************************
C
C    REMEMBER SOME OF THESE DIMENSIONS MUST ALSO BE CHANGED IN SUBROUTINES.
C
CHENN>
C      PARAMETER (NMAX=6000)
C      PARAMETER (ISPOTS=140)
C      PARAMETER (IRAST=80)
C      PARAMETER (MAXDIM=4096)
C      PARAMETER (ISLICE=8)
C
      PARAMETER (NMAX=48000)
      PARAMETER (ISPOTS=140)
      PARAMETER (IRAST=80)
      PARAMETER (MAXDIM=15000)
      PARAMETER (ISLICE=8)
C
      PARAMETER (LOOKUPMAX=10240)
CHENN<
      DIMENSION XCOORD(NMAX),YCOORD(NMAX),SINT(NMAX),BACK(NMAX),
     1ICORR(NMAX),XCOREF(NMAX),YCOREF(NMAX)
      INTEGER*2 NELH(8),NELK(8)
      INTEGER*4 JBLEFT(6)
      INTEGER*2 LJSCAN(6),JBACK(6),IPROP(NMAX),JBPASS(6),JBDONE(6)
      INTEGER*4  NX,NY
      INTEGER*2  JSTRIP(ISLICE,MAXDIM),JSTORE(ISPOTS,IRAST,IRAST)
      INTEGER*2 NBDBCK(NMAX),NGDBCK(NMAX)
      INTEGER*2 JH(NMAX),JK(NMAX),IFLAG(NMAX),IOVER(NMAX),ISTORE(NMAX)
      INTEGER*2 XBACK(6,NMAX),YBACK(6,NMAX)
      INTEGER*2 XBSTR(6),YBSTR(6)
      DIMENSION TOPT(10),BOTT(10),ITABLE(10,12),ITAV(12)
      DIMENSION NRAD(8),THETA(8),DELRAD(8),YCURVE(MAXDIM)
      DIMENSION TITLE(18),ODBACK(MAXDIM*2),STNDEV(MAXDIM*2)
CTSH      DIMENSION NXYZ(3),MXYZ(3),TEXT(20),IEXT(7),ARRAY(MAXDIM)
CTSH++
      DIMENSION NXYZ(3),MXYZ(3),IEXT(7),ARRAY(MAXDIM)
      CHARACTER*80 TEXT
CTSH--
      DIMENSION N9XYZ(3),M9XYZ(3),IEXT9(12)
CHENN>
C      REAL*4 LOOKUP(1500),LAMBDA
      REAL*4 LOOKUP(LOOKUPMAX),LAMBDA
CHENN<
      LOGICAL LPRINT
C
      PI=3.141592654
      MDONE=0
      NCYCLE =1
      ISTOP =0
CHENN>
C      IOVERL=700
C      IUNDER=1
C
      IOVERL=10000
      IUNDER=-30
CHENN<
C      ISIZE=2048
      MININT=-250
C      ISLICE=8
      WRITE(6,513)
C  READ IN CELL DIMENSIONS FOR USE BY AUTOMATIC TILT CALCULATION(EMTILT).
      READ(5,*)A,B,G,ANGDIS,SHRINK,TILTDIR,RKV,LPRINT
      XKV=RKV*1000.0
      LAMBDA=12.3/SQRT(XKV+XKV**2/(10.0**6.0))
      GMSTAR=180.0-G
      ASTAR=1.0/(A*SIN(GMSTAR*PI/180.0))
      BSTAR=1.0/(B*SIN(GMSTAR*PI/180.0))
      WRITE(6,9986)A,B,G,ANGDIS,SHRINK,TILTDIR,RKV,LAMBDA
9986  FORMAT(' CELL DIMENSIONS A,B,G',3F10.1/' ANGLE ON FILM FOR ',
     1'CORRECTION OF E.M. DISTORTION AND CORRECTION FACTOR',2F10.5/
     2' DIRECTION OF TILT ON MICROSCOPE (-1=CLOCKWISE) ',F6.1/
     3' Microscope voltage, wavelength.................',F6.1,F8.3)
      ANGDIS=ANGDIS*PI/180.0
C
      CALL IMOPEN(1,'IN','RO')
C
      CALL IRDHDR(1,NXYZ,MXYZ,MODE,DMIN,DMAX,DMEAN)
      NX=NXYZ(1)
      NY=NXYZ(2)
      IF(NX.GT.MAXDIM.OR.NY.GT.MAXDIM) THEN
      	WRITE(6,9982) NX,NY,MAXDIM
9982	FORMAT(' program dimensions too small - MAXDIM=',
     .		I5,'< NX,NY=',2I5)
      	STOP
      ENDIF
      ISIZE=NX
      SIZE=ISIZE
C         get extra parameters from actual image header
      CALL IRTEXT(1,IEXT,1,7)
      NSTEP=IEXT(1)/25.0
      WRITE(6,9985)(IEXT(J),J=1,7),NSTEP
9985  FORMAT(' IEXTRA(7) in PICKAUTO IN header',7I6,',  NSTEP=',I5)
      IF(NSTEP.EQ.0) NSTEP=4	! default for no header value
      STEP10=NSTEP*10
      X0 =IEXT(2)/STEP10 + NX/2.0
      Y0 =IEXT(3)/STEP10 + NY/2.0
      DX1=IEXT(4)/STEP10
      DY1=IEXT(5)/STEP10
      DX2=IEXT(6)/STEP10
      DY2=IEXT(7)/STEP10
C
      CALL IMOPEN(9,'INPARAM','RO')
      CALL IRDHDR(9,N9XYZ,M9XYZ,MODE9,D9MIN,D9MAX,D9MEAN)
      N9X=N9XYZ(1)
      N9Y=N9XYZ(2)
C         get extra parameters from header of autoindex image
      CALL IRTEXT(9,IEXT9,1,12)
      NSTEP=IEXT9(1)/25.0
      WRITE(6,9988)(IEXT9(J),J=1,12),NSTEP
9988  FORMAT(' IEXTRA(12) from INPARAM header ',12I6,',  NSTEP=',I5)
C         use autoindex parameters only if none in real image
      WRITE(6,9500) DX1,DY1,DX2,DY2
9500  FORMAT(' DX1,DY1,DX2,DY2 from input file',4F9.3)
      IF(SQRT(DX1**2+DY1**2).LT.1.0.OR.SQRT(DX2**2+DY2**2).LT.1.0)THEN
      	IF(NSTEP.EQ.0) NSTEP=4	! default for no header value
      	STEP10=NSTEP*10
      	DX1=IEXT9(4)/STEP10
      	DY1=IEXT9(5)/STEP10
      	DX2=IEXT9(6)/STEP10
      	DY2=IEXT9(7)/STEP10
      ENDIF
      WRITE(6,9501) DX1,DY1,DX2,DY2
9501  FORMAT(' DX1,DY1,DX2,DY2 from autoindex file',4F9.3)
C   
2002  FORMAT(' NX=',I5,'   NY=',I5)
      WRITE(6,2002) NX,NY
      WRITE(6,2001)NSTEP,DX1,DY1,DX2,DY2
2001  FORMAT(' INFORMATION FROM HEADER'/' NSTEP=',I5/' VECTOR (1,0)',
     12F10.2,';  VECTOR (0,1)',2F10.2/)
CHENN>
C      DO 1005 J=1,1500
C
      DO 1005 J=1,LOOKUPMAX
CHENN<
1005  LOOKUP(J)=0.
C     FIRST GENERATE LOOK-UP TABLE BY LINEAR INTERPOLATION
1003  READ(5,*) J,K
CHENN>
      WRITE(6,2992) J,K
2992  FORMAT(' J=',I5,'   K=',I5)
CHENN<
      LOOKUP(J)=K
      IF(J.EQ.1) GO TO 1003
      L=J
1006  L=L-1
      IF(LOOKUP(L).EQ.0.) GO TO 1006
      IF(L.EQ.(J-1)) GO TO 1008
      N=L+1
      DO 1007 M=N,J
1007  LOOKUP(M)=(LOOKUP(L)*(J-M)*2+J-L+LOOKUP(J)*(M-L)*2)/((J-L)*2)
CHENN>
C1008  IF(J.NE.1500) GO TO 1003
1008  IF(J.NE.LOOKUPMAX) GO TO 1003
CHENN<
      WRITE(6,1021)
1021  FORMAT('    ODMEAS    ODCORR----LOOK-UP TABLE TO CORRECT OD')
      DO 1009 J=1,1401,100
1009  WRITE(6,1020)J,LOOKUP(J)
1020  FORMAT(I10,F10.2)
1022  FORMAT(18A4)
1023  FORMAT('  LOOK-UP TABLE FOR ',18A4)
      READ(5,1022)TITLE
      WRITE(6,1023)TITLE
C    READ IN MEAN RADIAL BACKGROUND DENSITY
      READ(3,403) NCOMP,TITLE
      IMAX=1+SQRT(((NX+1.0)/2.0)**2+((NY+1.0)/2.0)**2)
C      WRITE(6,1024) IMAX
C1024  FORMAT(' IMAX =',I6)
      DO 1100 J=1,IMAX
      READ(3,509) ODBACK(J),STNDEV(J)
1100  CONTINUE
      WRITE(6,510) NCOMP,TITLE
      DO 1110 J=1,IMAX,100
      I=J-1
1110  WRITE(6,511) I,ODBACK(J),STNDEV(J)
C
C  READ IN MEAN Y-VARIATION CURVE FOR CORRECTION OF ALL DENSITY MEASUREMENTS.
      READ(4,1120)NYCOMP
      WRITE(6,1121)NYCOMP
1120  FORMAT(I10,F10.5)
1121  FORMAT(' Y-VARIATION CURVE READ IN, IDENTIFIER',I10,'  BRIEFLY')
      DO 1122 J=1,NY
      READ(4,1120)I,YCURVE(J)
      IF(I.EQ.J) GO TO 1122
      WRITE(6,1123)
1123  FORMAT(' Y-CURVE DATA IS WRONG?????????')
      STOP
1122  CONTINUE
      DO 1124 J=1,NY,200
1124  WRITE(6,1120)J,YCURVE(J)
      WRITE(6,399)
      READ(5,*) NPLATE
      READ(5,1022)TITLE
      WRITE(6,4049) NPLATE,TITLE
4049  FORMAT(' TITLE FOR OUTPUT',5X,I10,18A4)
      IF(NCOMP.NE.NPLATE) GO TO 590
      IF(NYCOMP.NE.NPLATE) GO TO 590
      READ(5,*) B1,B3
      WRITE(6,1002) B1,B3
1001  FORMAT(F10.5,F16.11)
1002  FORMAT(' STARTING VALUES FOR RADIAL DISTORTION PARAMETERS B1(NOT U
     1SED) AND B3  ARE',F10.5,F18.11)
      C1=0.0
      C3=0.0
      READ(5,*) CX0,CY0
      IF (CX0.NE.0.0.OR.CY0.NE.0.0) THEN
C     	  over-ride vectors on input diffraction pattern
      	X0=CX0
      	Y0=CY0
      	WRITE(6,8005) CX0,CY0
8005  	FORMAT(' Centre (X0,Y0) taken from card input',2F10.2)
      ENDIF
      READ(5,*) TLTAXA,TLTANG
      IF((TLTAXA.EQ.0.0).AND.(TLTANG.EQ.0.0)) GO TO 1025
      WRITE(6,1026)
1026  FORMAT(6(' ***'/),' *** WARNING, TLTAXA READ IN MUST BE ON FILM',
     16(' ***'/))
C
1025  READ(5,*)CDX1,CDY1,CDX2,CDY2
      IF (CDX1.EQ.0.0)GO TO 8002
C     OVER-RIDE VECTORS ON INPUT FILE
      DX1=CDX1
      DY1=CDY1
      DX2=CDX2
      DY2=CDY2
      WRITE(6,8001)
8001  FORMAT(' LATTICE VECTORS CORRECTED')
8002  WRITE(6,3)DX1,DY1,DX2,DY2
C
	IF(ABS(DX1+DX2+DY1+DY2).LE.0.0000001) THEN
	 WRITE(6,1029)
1029	 FORMAT(///' Lattice parameters not given either on input',
     .	 ' parameters or in header !!'//' !!!!!!!!!!!!!!!!!!!!!'///)
	 STOP
	ENDIF
C
      PI=3.141592654
      TLTAXA=TLTAXA*PI/180.0
      ANGA=ATAN2(DY1,DX1)
      ANGB=ATAN2(DY2,DX2)
      DANG=ANGB-ANGA
      IF(DANG.GT.PI) DANG=DANG-PI*2.0
      IF(DANG.LT.-PI) DANG=DANG+PI*2.0
      IF((TLTAXA.NE.0.0).OR.(TLTANG.NE.0.0)) GO TO 8003	! skip if tilt read in.
      ASTILT=SQRT(DX1**2+DY1**2)
      BSTILT=SQRT(DX2**2+DY2**2)
      GSTILT=DANG*180.0/PI
      IF(GSTILT.LT.0.0) GSTILT=-GSTILT
      CALL EMTILT(TL,TAXA,TANG,ASTAR,BSTAR,GMSTAR,ASTILT,BSTILT,GSTILT)
C TILTDIR DEFINED AT MICROSCOPE, -1 = CLOCKWISE ROTATION
C THESE CONVENTIONS APPLICABLE FOR 600 TO 1000 MM CAMERA LENGTH ON
C EM400 MICROSCOPE, WITH OBJECTIVE LENS CURRENT UNDERFOCUSSED
C TANG SHOULD ALWAYS BE RETURNED +VE FROM EMTILT
      TANG=TILTDIR*TANG
      WRITE(6,8004)TAXA,TANG,TL
8004  FORMAT(' TLTAXA AND TLTANG CALCULATED FROM LATTICE PARAMS AND',
     1' CELL DIMENSIONS IN PROGRAM     ',2F10.2,'   FILM ANGLE ',F10.2)
      TLTAXA=TL*PI/180.0
      TLTANG=TANG
8003  CONTINUE
      IF(DANG.LT.0.) THEN
      	TLTDIR=ANGA+TLTAXA
      ELSE
      	TLTDIR=ANGA-TLTAXA
      ENDIF
      IF(TLTDIR.GT.PI) TLTDIR=TLTDIR-PI
      IF(TLTDIR.LT.-PI) TLTDIR=TLTDIR+PI
      IF(TLTDIR.LT.0.) TLTDIR=TLTDIR+PI
C
C  THIS ENSURES TLTDIR IS BETWEEN 0 AND PI.
C  TLTNRM IS THEREFORE BETWEEN PI/2 AND 3*PI/2.
C  THUS, FOR CURRENT (JUNE 82) CONVENTION FOR FILM SCANNING AND ELECTRON
C  DIFFRACTION  TLTANG IS ALWAYS NEGATIVE.
C
      TLTNRM=TLTDIR+PI/2.0
      TLTNRM=TLTNRM*180.0/PI
      TLTDIR=TLTDIR*180.0/PI
      TLTAXA=TLTAXA*180.0/PI
      WRITE(6,1027)TLTDIR
1027  FORMAT(' ABSOLUTE DIRECTION OF TILTAXIS ON FILM (RELATIVE TO',
     $' X SCAN AXIS) WAS',19('-'),F7.1,'  DEGREES')
      TLTDIR=TLTDIR*PI/180.0
1     FORMAT(6F10.2)
      WRITE(6,2) X0,Y0,TLTAXA,TLTANG
      TLTNRM=TLTNRM*PI/180.0
      TLTANG=TLTANG*PI/180.0
      TLTAXA=TLTAXA*PI/180.0
2     FORMAT(' STARTING COORDS X0,Y0  ',62('-'),2F10.2/
     1' ANGLE BETWEEN ASTAR AND TILTAXIS ON FILM AND TLTANG',33('-'),
     2 2F10.2)
3     FORMAT(' INITIAL POSITION OF (1,0)  ',60('-'),2F10.3/
     1' INITIAL POSITION OF (0,1)  ',60('-'),2F10.3)
C
C     MAKE SCAN UP TO INDEX IN EACH HEXTANT OF 18 - - - 3.0 ANGSTROMS
C     BUT PUT IN A RESOLUTION LIMIT TO AVOID TOO MUCH DATA.
C     THIS MEANS WE NEED FROM -20 TO +20 IN HEXAGONAL LATTICE .
C
      READ(5,*) ROUT,RIN,PRPMAX	       ! now read radii in, in Angstroms
		WRITE(6,1033) RIN,ROUT,PRPMAX
1033		FORMAT(' Inner,outer(para/perp) radii input',30('-'),3F6.1)
C	convert to radius in pixels, assuming ASTAR corresponds to (1,0)
		AH     = SQRT(DX1**2+DY1**2)
		AHCORR = (COS(TLTAXA))**2+(SIN(TLTAXA)*COS(TLTANG))**2
		AHCORR = AH*SQRT(AHCORR)
		ROUT   = AHCORR/ASTAR*(1.0/ROUT)
		RIN    = AHCORR/ASTAR*(1.0/RIN)
		PRPMAX = AHCORR/ASTAR*(1.0/PRPMAX)
      		ROUTSTATS=ROUT
C	now in pixels
C
      READ(5,*) NXM,NYM,NSM,NXB,NYB,NSB,NPOS,
     1NXMT,NYMT,NSMT,NXBT,NYBT,NSBT,NPOST,NTYPE
      IF(NXM.EQ.0) THEN
C          over-ride raster sizes on input files
		NXM=IEXT9(8)
		NYM=IEXT9(9)
		NXMT=0
		NYMT=0
		NXB=IEXT9(10)
		NYB=IEXT9(11)
		NTYPE=IEXT9(12)
		NSM=1
		NSB=1
		NPOS=1
      ENDIF
      IF((NXM.GT.IRAST).OR.(NYM.GT.IRAST)) GO TO 1032
      IF((NXMT.LE.IRAST).AND.(NYMT.LE.IRAST)) GO TO 1030
1032  WRITE(6,1031)NXM,NYM,NXMT,NYMT
      STOP
1031  FORMAT(' MEASURING RASTER TOO BIG FOR STORE',4I5)
1030  IF(NPOS.EQ.0) NPOS=1
      IF(NPOST.EQ.0) NPOST=1
6     FORMAT(15I5)
      READ(5,*) FRACT,ABSOL,XAMINE,NCYC,NCYC1
C
      READ(5,*)NELIM
      IF(NELIM.EQ.0)GO TO 35
      IF(NELIM.GT.0)GO TO 1200
      MELIM=-NELIM
      READ(5,*)(NELH(I),NELK(I),I=1,MELIM)
      WRITE(6,1201)(NELH(I),NELK(I),I=1,MELIM)
1201  FORMAT(' REFLECTIONS WHOSE RASTERS ARE TO BE PRINTED FOR',
     1' INSPECTION'/(2X,2I4))
      GO TO 35
1200  READ(5,*)(NELH(I),NELK(I),I=1,NELIM)
      WRITE(6,36)(NELH(I),NELK(I),I=1,NELIM)
36    FORMAT(' REFLECTIONS TO BE EXCLUDED FROM LATTICE REFINEMENT'/
     1(2X,2I4))
35    IF(NCYC.EQ.0)ISTOP=1
C
C
C    STORAGE LIMITS RASTER SIZE TO IRAST X IRAST.
C    MEASURING AND BACKGROUND RASTERS - - MAKE BACKGROUND AS BIG AS POSSIBLE
C
      WRITE(6,4) ROUT,RIN,PRPMAX
4     FORMAT(' DATA INCLUDED OUT TO RADIUS',57('-'),F10.2,'   FROM',
     1F10.2/' WITH MAXIMUM PERPENDICULAR DISTANCE FROM TILT AXIS',
     2 35('-'),F10.2)
      WRITE(6,15) NXM,NYM,NSM,NXB,NYB,NSB,NPOS
15    FORMAT(' MEASURING AND BACKGROUND RASTERS AND POSITIONS ARE',
     1I3,' X',I3,' X',I3,' AND',I3,' X',I3,' X',I3,' NPOS=',I1)
      IF(NXMT.NE.0) GO TO 30
      NXMT=NXM
      NYMT=NYM
      NSMT=NSM
      NXBT=NXB
      NYBT=NYB
      NSBT=NSB
      NPOST=NPOS
30    CONTINUE
      WRITE(6,31) NXMT,NYMT,NSMT,NXBT,NYBT,NSBT,NPOST
31    FORMAT(' AND FAR AWAY FROM TILT AXIS ARE',17X,
     1I3,' X',I3,' X',I3,' AND',I3,' X',I3,' X',I3,' NPOS=',I1)
C    BACKGROUND POSITION SHOULD BE NEARER SPOTS WHEN THE BACKGROUND
C    DENSITY IS CURVED. TRY NPOS=1 FOR BIG SPOTS, NPOS=4 FOR SMALL.
      WRITE(6,32)NTYPE
32    FORMAT(' TYPE OF BACKGROUND POSITIONS',20X,I3)
      IF(NTYPE.EQ.0.OR.NTYPE.EQ.1)NUMIB=6
      IF(NTYPE.EQ.2)NUMIB=4
C
C
C     GENERATE AND SORT COORDINATES OF ALL SPOTS OUT TO REQUIRED RESOLUTION.
260   NDATA=0
      IF(ISTOP.EQ.0) THEN
        IF(LPRINT) WRITE(6,9028)
9028    FORMAT(//' NEW CYCLE'//)
      ENDIF
      IF(ISTOP.EQ.1) THEN
        IF(LPRINT) WRITE(6,9060)
9060  FORMAT(//' FINAL COMPUTATION OF BACKGROUND-CORRECTED INTENSITIES'
     1//)
      ENDIF
C
      AH=SQRT(DX1**2+DY1**2)
      AHCORR=(COS(TLTAXA))**2+(SIN(TLTAXA)*COS(TLTANG))**2
      AHCORR=AH*SQRT(AHCORR)
C  only one astar/ahcorr needed because EWALD is required to be in pixels
      EWALD=(LAMBDA*0.5*TAN(TLTANG)*ASTAR/AHCORR)
      XEWALD=EWALD*COS(TLTNRM)
      YEWALD=EWALD*SIN(TLTNRM)
      CORMAX=EWALD*ROUT*ROUT
      WRITE(6,263) CORMAX,ROUT
263   FORMAT(' MAXIMUM EWALD SPHERE POSITION CORRECTION=',F7.2,
     1'   AT RADIUS',F9.2)
C 
      IF(NTYPE.NE.0)GO TO 33
C     NTYPE=0
      XBA=(DX1+DX2)/3
      YBA=(DY1+DY2)/3
      XBB=(2*DX1-DX2)/3
      YBB=(2*DY1-DY2)/3
      XBC=(DX1-2*DX2)/3
      YBC=(DY1-2*DY2)/3
      GO TO 34
33    IF (NTYPE.NE.1)GO TO 37
C     NTYPE=1
      XBA=(2*DX1+DX2)/3
      YBA=(2*DY1+DY2)/3
      XBB=(DX1+2*DX2)/3
      YBB=(DY1+2*DY2)/3
      XBC=(DX1-DX2)/3
      YBC=(DY1-DY2)/3
      GO TO 34
C     NTYPE=2; 4 BACKGROUNDS 
37    XBA=(DX1+DX2)/2
      YBA=(DY1+DY2)/2
      XBB=(DX1-DX2)/2
      YBB=(DY1-DY2)/2
      XBC=0.
      YBC=0.
34    CONTINUE
C
C
      DO 261 K=1,ISPOTS
261   JSTORE(K,1,1)=-1
      DO 20 I=1,ISPOTS+1
      DO 20 J=1,ISPOTS+1
      IH=I-1-ISPOTS/2
      IK=J-1-ISPOTS/2
      RSQ=(IH*DX1+IK*DX2)**2+(IH*DY1+IK*DY2)**2
      XCSTRT=X0+IH*DX1+IK*DX2+RSQ*XEWALD
      YCSTRT=Y0+IH*DY1+IK*DY2+RSQ*YEWALD
      XC=XCSTRT
      YC=YCSTRT
      DO 262 NAP=1,3
      RADSQ=(XC-X0)**2+(YC-Y0)**2
      XB3CR=(XC-X0)*B3*RADSQ
      YB3CR=(YC-Y0)*B3*RADSQ
      XC=XCSTRT+XB3CR
      YC=YCSTRT+YB3CR
262   CONTINUE
      DXC=XC-X0
      DYC=YC-Y0
      RAD=SQRT(DXC*DXC+DYC*DYC)
      IF((RAD.GT.ROUT).OR.(RAD.LT.RIN)) GO TO 20
      IF(XC.GT.FLOAT(NX)) GO TO 20
      IF(XC.LT.1.0) GO TO 20
      IF(YC.GT.FLOAT(NY)) GO TO 20
      IF(YC.LT.1.0) GO TO 20
10    CONTINUE
      L=0
C    CALCULATION OF DISTANCE FROM TILTAXIS AND APPROPRIATE
C    SIZE OF SPOT AND BACKGROUND RASTER SIZES AND RASTER POSITIONS.
      DPERP=DYC*COS(TLTDIR)-DXC*SIN(TLTDIR)
      IF(DPERP.LT.0.) DPERP=-DPERP
      IF (DPERP.GT.PRPMAX)GO TO 20
C
      PROPOR=DPERP/PRPMAX
      NXSH=0.5*(PROPOR*(NXMT-NXM)+NXM)
      NYSH=0.5*(PROPOR*(NYMT-NYM)+NYM)
      IF(DXC.LT.0)NXSH=-NXSH
      IF(DYC.LT.0)NYSH=-NYSH
      REXTR=(DXC*NXSH+DYC*NYSH)/RAD
      IF((RAD+REXTR).GT.ROUT)GO TO 20
      DIV=1.0/(PROPOR*(NPOST-NPOS)+NPOS)
      DIVDR=DIV/RAD
      NXBTH=0.5*((NXBT-NXB)*PROPOR+NXB)/DIV
      NYBTH=0.5*((NYBT-NYB)*PROPOR+NYB)/DIV
      IF(XBA.GE.0)XBAT=XBA+NXBTH
      IF(XBA.LT.0)XBAT=XBA-NXBTH
      IF(XBB.GE.0)XBBT=XBB+NXBTH
      IF(XBB.LT.0)XBBT=XBB-NXBTH
      IF(XBC.GE.0)XBCT=XBC+NXBTH
      IF(XBC.LT.0)XBCT=XBC-NXBTH
      IF(YBA.GE.0)YBAT=YBA+NYBTH
      IF(YBA.LT.0)YBAT=YBA-NYBTH
      IF(YBB.GE.0)YBBT=YBB+NYBTH
      IF(YBB.LT.0)YBBT=YBB-NYBTH
      IF(YBC.GE.0)YBCT=YBC+NYBTH
      IF(YBC.LT.0)YBCT=YBC-NYBTH
      REXT1=(DXC*XBAT+DYC*YBAT)*DIVDR
      REXT2=(DXC*XBBT+DYC*YBBT)*DIVDR
      REXT3=(DXC*XBCT+DYC*YBCT)*DIVDR
      IF((RAD+REXT1).GT.ROUT)L=L+1
      IF((RAD-REXT1).GT.ROUT)L=L+2
      IF((RAD+REXT2).GT.ROUT)L=L+4
      IF((RAD-REXT2).GT.ROUT)L=L+8
      IF((RAD+REXT3).GT.ROUT.AND.NTYPE.NE.2)L=L+16
      IF((RAD-REXT3).GT.ROUT.AND.NTYPE.NE.2)L=L+32
      IF(NTYPE.EQ.2)L=L+48
      NDATA=NDATA+1
      IF(NDATA.LE.NMAX)GO TO 9010
      WRITE(6,7)
7     FORMAT('  *****ERROR - TOO MANY SPOTS FOR STORE')
      STOP
9010  CONTINUE
      IPROP(NDATA)=PROPOR*1000.
      JH(NDATA)=IH
      JK(NDATA)=IK
      XCOORD(NDATA)=XC
      YCOORD(NDATA)=YC
      SINT(NDATA)=0.0
      BACK(NDATA)=0.0
      NGDBCK(NDATA)=0
      NBDBCK(NDATA)=0
      IFLAG(NDATA)=0
      IOVER(NDATA)=0
      XCOREF(NDATA)=0.0
      YCOREF(NDATA)=0.0
      ICORR(NDATA)=0.0
      ISTORE(NDATA)=-9999
      XBACK(1,NDATA)=XC+XBA*DIV
      XBACK(2,NDATA)=XC-XBA*DIV
      XBACK(3,NDATA)=XC+XBB*DIV
      XBACK(4,NDATA)=XC-XBB*DIV
      XBACK(5,NDATA)=XC+XBC*DIV
      XBACK(6,NDATA)=XC-XBC*DIV
      YBACK(1,NDATA)=YC+YBA*DIV
      YBACK(2,NDATA)=YC-YBA*DIV
      YBACK(3,NDATA)=YC+YBB*DIV
      YBACK(4,NDATA)=YC-YBB*DIV
      YBACK(5,NDATA)=YC+YBC*DIV
      YBACK(6,NDATA)=YC-YBC*DIV
C   SET BACKGROUND CORRDS TO ZERO IF OUTSIDE RANGE
      IF(L.EQ.0)GO TO 20
      IF(L.NE.1.AND.L.NE.3.AND.L.NE.5.AND.L.NE.9.AND.L.NE.17.AND.L.NE.33
     1)GO TO 9021
      XBACK(1,NDATA)=0
      YBACK(1,NDATA)=0
      L=L-1
      IF(L.EQ.0)GO TO 20
9021  IF(L.NE.2.AND.L.NE.6.AND.L.NE.10.AND.L.NE.18.AND.L.NE.34)GO TO 902
     12
      XBACK(2,NDATA)=0
      YBACK(2,NDATA)=0
      L=L-2
      IF(L.EQ.0)GO TO 20
9022   IF(L.NE.4.AND.L.NE.12.AND.L.NE.20.AND.L.NE.36)GO TO 9023
      XBACK(3,NDATA)=0
      YBACK(3,NDATA)=0
      L=L-4
      IF(L.EQ.0)GO TO 20
9023  IF(L.NE.8.AND.L.NE.24.AND.L.NE.40)GO TO 9024
      XBACK(4,NDATA)=0
      YBACK(4,NDATA)=0
      L=L-8
      IF (L.EQ.0)GO TO 20
9024  IF(L.NE.16.AND.L.NE.48)GO TO 9025
      XBACK(5,NDATA)=0
      YBACK(5,NDATA)=0
C     WRITE(6,10005)IH,IK,XC,YC
C10005 FORMAT(' IB=5 IS OUT FOR',2I5,2F10.2)
      L=L-16
      IF(L.EQ.0)GO TO 20
9025  IF(L.NE.32)GO TO 9026
      XBACK(6,NDATA)=0
      YBACK(6,NDATA)=0
C     WRITE(6,10006)IH,IK,XC,YC
C10006 FORMAT(' IB=6 IS OUT FOR',2I5,2F10.2)
C
C
      GO TO 20
C  MORE THAN TWO BACKGROUNDS OUTSIDE LIMIT
9026   NDATA=NDATA-1
C     WRITE(6,9027)IH,IK,DXC,DYC,L,RAD,REXT1,REXT2,REXT3
C9027  FORMAT(' BACKGROUNDS OUTSIDE MAX RADIUS',2I5,2F8.2,I5,4F8.2)
C
20    CONTINUE
      WRITE(6,5) NDATA
5     FORMAT(//' NUMBER OF SPOTS WITHIN LIMITS = ',I10)
C     SORT SPOTS IN ORDER OF INCREASING Y-COORDINATE
      NDATAB = NDATA-1
      DO 50 I=1,NDATAB
      JST = I+1
      DO 50 J=JST,NDATA
      IF(YCOORD(I).LE.YCOORD(J)) GO TO 50
      YSTORE=YCOORD(I)
      XSTORE=XCOORD(I)
      JHSTOR=JH(I)
      JKSTOR=JK(I)
      IPRPST=IPROP(I)
      YCOORD(I)=YCOORD(J)
      XCOORD(I)=XCOORD(J)
      JH(I)=JH(J)
      JK(I)=JK(J)
      IPROP(I)=IPROP(J)
      YCOORD(J)=YSTORE
      XCOORD(J)=XSTORE
      JH(J)=JHSTOR
      JK(J)=JKSTOR
      IPROP(J)=IPRPST
C
      DO 9050 K=1,NUMIB
      XBSTR(K)=XBACK(K,I)
      YBSTR(K)=YBACK(K,I)
      XBACK(K,I)=XBACK(K,J)
      YBACK(K,I)=YBACK(K,J)
      XBACK(K,J)=XBSTR(K)
      YBACK(K,J)=YBSTR(K)
9050  CONTINUE
C
50    CONTINUE
C
      JSUM=(NDATA*(NDATA+1))/2
C
C
C
C     READ IN ISLICE LINES AT A TIME AND BUILD UP A RASTER ROUND EACH SPOT
C     TO BE MEASURED .
C     USED TO BE 16 LINES, SO VARIABLE IS LINE16.
C
      LINE1=1
      LINE16=ISLICE
      JBEGIN=1
      JPASS=0
      JSKIP=0
      JDONE=0
      JLEFT=JSUM
      LSCAN=0
      JSTART=1
      DO 61 N=1,NUMIB
      LJSCAN(N)=0
      JBPASS(N)=0
      JBDONE(N)=0
      JBLEFT(N)=JSUM
61    JBACK(N)=1
65    DO 70 L=1,ISLICE
      CALL IRDLIN(1,ARRAY,*8999)
      DO 67 M=1,NX
67    JSTRIP(L,M)=ARRAY(M)
70    CONTINUE
      IF(JSTART.GT.NDATA) GO TO 81
      IOK=1
      JCHECK=JSTART-ISPOTS
      IF(JCHECK.LT.1) JCHECK=1
      DO 72 JCH=JCHECK,JSTART
      IF(ISTORE(JCH).LT.1) GO TO 72
      JCHECK=JCH
      GO TO 73
72    CONTINUE
73    JSCAN=JSTART+ISPOTS
      IF(JSCAN.GT.NDATA) JSCAN=NDATA
      DO 80 J=JCHECK,JSCAN
      IF(J.EQ.NDATA) IOK=0
      IF(ISTORE(J).EQ.-9000) GO TO 80
      PROP=IPROP(J)/1000.
      NXMRA=PROP*(NXMT-NXM)+NXM
      NYMRA=PROP*(NYMT-NYM)+NYM
      XMHALF=(NXMRA-1.0)/2.0
      YMHALF=(NYMRA-1.0)/2.0
      NBEGIN=YCOORD(J)-YMHALF+0.5
      NFINSH=NBEGIN+NYMRA-1
      NXST=XCOORD(J)-XMHALF+0.5
      NXFI=NXST+NXMRA-1
      IF((NXST.GE.1).AND.(NXFI.LE.NX)) GO TO 974
      IF(J.LE.LSCAN)GO TO 80
      IF(LPRINT) WRITE(6,5005)J,JH(J),JK(J)
5005  FORMAT(' X RASTER GOES OUTSIDE RANGE',3I5)
      ISTORE(J)=-9000
      JDONE=JDONE+1
      JLEFT=JLEFT-J
      JSKIP=JSKIP+1
      JPASS=J
      GO TO 80
974   IF(NBEGIN.GE.1) GO TO 75
      WRITE (6,5003)J,JH(J),JK(J)
5003  FORMAT(' NBEGIN <1 ',3I5)
      JSTART=J+1
      ISTORE(J)=-9000
      JDONE=JDONE+1
      JLEFT=JLEFT-J
      JSKIP=JSKIP+1
      JPASS=J
      GO TO 80
75    IF(NFINSH.LE.MAXDIM)GO TO 975
      ISTORE(J)=-9000
      IF(LPRINT) WRITE(6,5004)J,JH(J),JK(J)
5004  FORMAT(' NFINSH >MAXDIM ',3I5)
      JDONE=JDONE+1
      JLEFT=JLEFT-J
      JSKIP=JSKIP+1
      JPASS=J
      GO TO 80
975   IF(NBEGIN.GE.LINE1) GO TO 76
      IF(NFINSH.GE.LINE1) GO TO 77
      IF(LPRINT) WRITE(6,777)
      GO TO 80
77    IS=ISTORE(J)
      NFIN=NFINSH
      IF(NFINSH.GT.LINE16) NFIN=LINE16
C
       DO 85 IY=LINE1,NFIN
        DO 85 IX=NXST,NXFI
          IXRAST=IX+1-NXST
          IYRAST=IY+1-NBEGIN 
          IL=IY+1-LINE1 
          JSTORE(IS,IXRAST,IYRAST)=JSTRIP(IL,IX)
          JLOOK=JSTRIP(IL,IX)
CHENN>
C         IF(JLOOK.LT.1024.AND.JLOOK.GT.0)GO TO 6999
C           IF(LPRINT) WRITE(6,7000)JH(J),JK(J)
C7000       FORMAT(' NONSENSE OD  IN REFLECTION',2I5)
C
          IF(JLOOK.LT.LOOKUPMAX.AND.JLOOK.GT.0)GO TO 6999
            IF(LPRINT) WRITE(6,7000)JH(J),JK(J),JLOOK
7000        FORMAT(' NONSENSE OD  IN REFLECTION',3I5)
CHENN<
            IF(JLOOK.GE.1024)IOVER(J)=1
            IF(JLOOK.LE.1)IOVER(J)=2
          GO TO 85
6999      IF(JLOOK.GT.IOVERL) IOVER(J)=1
      IF(JLOOK.LE.IUNDER) IOVER(J)=2
      SINT(J)=SINT(J)+LOOKUP(JLOOK)
85    CONTINUE
C
      IF(NFINSH.GT.LINE16) GO TO 80
      JSTART=J+1
      JDONE=JDONE+1
      JLEFT=JLEFT-J
      IPRINT =0
C     WRITE(6,5001) J,JH(J),JK(J),IPROP(J),ISTORE(J),NXMRA,LINE1
      IF((J/20)*20.EQ.J) IPRINT=1
      IF(NCYCLE.EQ.1.AND.NCYC1.NE.0) IPRINT=1
      IF(.NOT.LPRINT) IPRINT=0
C
      IF(NELIM.EQ.0)GO TO 5002
      IF(NELIM.GT.0)GO TO 10012
C     PRINT SELECTED RASTERS
      DO 9231 IELIM=1,MELIM
      IF(JH(J).NE.NELH(IELIM).OR.JK(J).NE.NELK(IELIM))GO TO 9231
      IF(LPRINT) THEN
      	IPRINT=2
      ELSE
      	IPRINT=0
      ENDIF
      GO TO 5002
9231  CONTINUE
10012 DO 5000 IELIM=1,NELIM
      IF(JH(J).NE.NELH(IELIM).OR.JK(J).NE.NELK(IELIM))GO TO 5000
      ABNORM=-1.0
      IF(LPRINT) WRITE(6,5001)JH(J),JK(J)
5001  FORMAT(' REFLECTION EXCLUDED FROM LATTICE REFINEMENT',2I5)
      GO TO 5002
5000  CONTINUE
C
5002  CALL CENTRE(JSTORE,IS,XCEN,YCEN,ABNORM,IPRINT,NXMRA,NYMRA,
     1JH(J),JK(J),XAMINE,ODBACK,ODBTOT,NXST,NBEGIN,X0,Y0,
     2ROUT,RIN,YCURVE,
     3LOOKUP,IOVERL,IUNDER,IOVERC,ISTOP)
      ISTORE(J)=-9000
      IF(ISTOP.EQ.1) then
        TCORR=SINT(J)-ODBTOT
        ICORR(J)=NINT(TCORR)
        IF(IOVERC.EQ.1)IOVER(J)=1
        IF(IOVERC.EQ.2)IOVER(J)=2
        IF(IOVERC.EQ.3)IOVER(J)=3
        IF(ABNORM.NE.0.0)GO TO 80
        XCOREF(J)=XCEN+NXST-1.0
        YCOREF(J)=YCEN+NBEGIN-1.0
        GO TO 80
      end if
      IF(ABNORM.NE.0.0) GO TO 80
      XCOREF(J)=XCEN+NXST-1.0
      YCOREF(J)=YCEN+NBEGIN-1.0
      GO TO 80
76    IF(NBEGIN.LE.LINE16) GO TO 74
      IOK=0
      GO TO 80
74    DO 78 K=1,ISPOTS
      IF(JSTORE(K,1,1).EQ.-1) GO TO 79
78    CONTINUE
      WRITE(6,82)ISPOTS
82    FORMAT(//' RASTER STORAGE TOO SMALL, NEEDS .GT.',I5)
      STOP
79    IS=K
      ISTORE(J)=K
      NFIN=NFINSH
      IF(NFINSH.GT.LINE16) GO TO 984
      GO TO 988
984   NFIN=LINE16
988   DO 83 IY=NBEGIN,NFIN
      DO 83 IX=NXST,NXFI
      IXRAST=IX+1-NXST
      IYRAST=IY+1-NBEGIN
      IL=IY+1-LINE1
      JSTORE(IS,IXRAST,IYRAST)=JSTRIP(IL,IX)
      JLOOK=JSTRIP(IL,IX)
CHENN>
C      IF(JLOOK.LT.1024.AND.JLOOK.GT.0)GO TO 6998
C      IF(LPRINT) WRITE(6,7000)JH(J),JK(J)
C
      IF(JLOOK.LT.LOOKUPMAX.AND.JLOOK.GT.0)GO TO 6998
      IF(LPRINT) WRITE(6,7000)JH(J),JK(J),JLOOK
CHENN<
      IF(JLOOK.GE.1024)IOVER(J)=1
      IF(JLOOK.LE.1)IOVER(J)=2
      GO TO 83
6998  IF(JLOOK.GT.IOVERL) IOVER(J)=1
      IF(JLOOK.LE.IUNDER) IOVER(J)=2
      SINT(J)=SINT(J)+LOOKUP(JLOOK)
 83   CONTINUE
      IF(NFINSH.GT.LINE16) GO TO 80
      JSTART=J+1
      JDONE=JDONE+1
      JLEFT=JLEFT-J
C     PRINT OUT EVERY 20TH RASTER
      IPRINT=0
C     WRITE(6,5002) J,JH(J),JK(J),IPROP(J),ISTORE(J),NXMRA,LINE1
      IF((J/20)*20.EQ.J) IPRINT=1
      IF(NCYCLE.EQ.1.AND.NCYC1.NE.0) IPRINT=1
      IF(.NOT.LPRINT) IPRINT=0
C
      IF (NELIM.EQ.0)GO TO 10002
      IF(NELIM.GT.0)GO TO 10010
C     PRINT SELECTED RASTERS
      DO 10011 IELIM=1,MELIM
      IF(JH(J).NE.NELH(IELIM).OR.JK(J).NE.NELK(IELIM))GO TO 10011
      IF(LPRINT) THEN
      	IPRINT=2
      ELSE
      	IPRINT=0
      ENDIF
      GO TO 10002
10011 CONTINUE
10010 DO 10001 IELIM=1,NELIM
      IF(JH(J).NE.NELH(IELIM).OR.JK(J).NE.NELK(IELIM))GO TO 
     1 10001
      ABNORM=-1.0
      IF(LPRINT) WRITE(6,5001)JH(J),JK(J)
      GO TO 10002
10001 CONTINUE
C      
10002 CALL CENTRE(JSTORE,IS,XCEN,YCEN,ABNORM,IPRINT,NXMRA,NYMRA,
     1JH(J),JK(J),XAMINE,ODBACK,ODBTOT,NXST,NBEGIN,X0,Y0,
     2ROUT,RIN,YCURVE,
     3LOOKUP,IOVERL,IUNDER,IOVERC,ISTOP)
      ISTORE(J)=-9000
      IF(ISTOP.EQ.1)  THEN
  	TCORR=SINT(J)-ODBTOT
      	ICORR(J)=NINT(TCORR)
      	IF(IOVERC.EQ.1)IOVER(J)=1
      	IF(IOVERC.EQ.2)IOVER(J)=2
      	IF(IOVERC.EQ.3)IOVER(J)=3
      	IF(ABNORM.NE.0.0)GO TO 80
      	XCOREF(J)=XCEN+NXST-1.0
      	YCOREF(J)=YCEN+NBEGIN-1.0
      	GO TO 80
      ENDIF
      IF(ABNORM.NE.0.0) GO TO 80
      XCOREF(J)=XCEN+NXST-1.0
      YCOREF(J)=YCEN+NBEGIN-1.0
80    CONTINUE
      LSCAN=JSCAN
      IF(IOK.EQ.0) GO TO 81
      WRITE(6,82) ISPOTS
      WRITE(6,82) ISPOTS
      STOP
81    CONTINUE
C
C     NOW DO EACH OF THE SIX BACKGROUND MEASUREMENTS .
C     SAME AS FOR MEASURING RASTER BUT NO STORAGE OF RASTER NUMBERS .
C
C     BACKGROUND RASTER INTEGRATION
C
C     SET IFLAG(J)=N   IF N BACKGROUND RASTERS ARE MISSED
C
      IF(ISTOP.EQ.0) GO TO 222
      DO 201 IB=1,NUMIB
      JINIT=JBACK(IB)
      IOK=1
      JCHECK=JINIT-ISPOTS
      JSCAN=JINIT+ISPOTS
      IF(JCHECK.LT.1) JCHECK=1
      IF(JSCAN.GT.NDATA) JSCAN=NDATA
      DO 200 J=JCHECK,JSCAN
      IF(J.EQ.NDATA) IOK=0
      PROP=IPROP(J)/1000.
      NXBRA=PROP*(NXBT-NXB)+NXB
      NYBRA=PROP*(NYBT-NYB)+NYB
      XBHALF=(NXBRA-1.0)/2.0
      YBHALF=(NYBRA-1.0)/2.0
      NBEGIN=YBACK(IB,J)-YBHALF +0.5
      NFINSH=NBEGIN + NYBRA-1
      NXST = XBACK(IB,J) - XBHALF + 0.5
      NXFI = NXST + NXBRA -1
C    SET FLAG IF UNSUITABLE BACKGROUND POSITION
      IF(XBACK(IB,J).EQ.0.AND.YBACK(IB,J).EQ.0)GO TO 101
      IF((NXST.LT.1).OR.(NXFI.GT.NX)) GO TO 101
      IF(NBEGIN.GE.1) GO TO 175
C     SET FLAG ON FIRST SLICE ONLY
      IF(LINE1.NE.1)GO TO 200
      IFLAG(J)=IFLAG(J)+1
      JTEST=JBPASS(IB)
      JBPASS(IB)=MAX0(J,JTEST)
      JBDONE(IB)=JBDONE(IB)+1
      JBLEFT(IB)=JBLEFT(IB)-J
      GO TO 200
C
C     FLAG UNLESS ALREADY FLAGGED ON PREVIOUS SLICE
101   IF(J.LE.LJSCAN(IB))GO TO 200
      IFLAG(J)=IFLAG(J)+1
      JTEST=JBPASS(IB)
      JBPASS(IB)=MAX0(J,JTEST)
      JBDONE(IB)=JBDONE(IB)+1
      JBLEFT(IB)=JBLEFT(IB)-J
      GO TO 200
175   CONTINUE
      IF(NFINSH.LE.MAXDIM)GO TO 9175
      IF(J.LE.LJSCAN(IB))GO TO 200
      IFLAG(J)=IFLAG(J)+1
      JTEST=JBPASS(IB)
      JBPASS(IB)=MAX0(J,JTEST)
C     WRITE(6,10012)IB,JH(J),JK(J),XBACK(IB,J),YBACK(IB,J),
C    1IB,JBPASS(IB)
      JBDONE(IB)=JBDONE(IB)+1
      JBLEFT(IB)=JBLEFT(IB)-J
      GO TO 200
9175  CONTINUE
      IF(NBEGIN.GE.LINE1) GO TO 176
      IF(NFINSH.GE.LINE1) GO TO 177
      IF(J.LT.JINIT) GO TO 200
      IF(LPRINT) WRITE(6,778) IB,J,JH(J),JK(J),IFLAG(J),
     $NBEGIN,NFINSH,LINE1,XBACK(IB,J),YBACK(IB,J)
778   FORMAT(' FUNNY ONE',8I5,2I10)
      GO TO 200
177   NFIN=NFINSH
      IF(NFINSH.GT.LINE16) NFIN=LINE16
      DO 185 IY=LINE1,NFIN
      IL=IY+1-LINE1
      YCO=IY-Y0
      YCOSQ=YCO**2
      DO 185 IX=NXST,NXFI
      JLOOK=JSTRIP(IL,IX)
CHENN>
C       IF(JLOOK.LT.1024.AND.JLOOK.GT.0)GO TO 6997
C       IF(LPRINT) WRITE(6,7002)JH(J),JK(J)
C 7002  FORMAT(' NONSENSE OD  IN BACKGROUND OF REFLECTION',2I5)
C
      IF(JLOOK.LT.LOOKUPMAX.AND.JLOOK.GT.0)GO TO 6997
      IF(LPRINT) WRITE(6,7002)JH(J),JK(J),JLOOK
7002  FORMAT(' NONSENSE OD  IN BACKGROUND OF REFLECTION',3I5)
CHENN<
      IF(JLOOK.GE.1024)IOVER(J)=1
      IF(JLOOK.LE.1)IOVER(J)=2
      GO TO 185
6997  IF(JLOOK.GT.IOVERL) IOVER(J)=1
      IF(JLOOK.LE.IUNDER) IOVER(J)=2
C  CORRECT BACKGROUND RASTER ODS WITH RADIAL DISTRIBUTION
      XCO=IX-X0
      RAD=1.+SQRT(YCOSQ+XCO**2)
      IRAD=RAD
      IF(RAD.LE.RIN)IOVER(J)=3
      IF (RAD.GT.(ROUT+1.))IOVER(J)=3
      DRAD=RAD-IRAD
      ODB=ODBACK(IRAD)*(1-DRAD)+ODBACK(IRAD+1)*DRAD+YCURVE(IY)
C     CHECK THAT OD IN BACKGROUND IS NOT TOO FAR AWAY FROM MEAN VALUE
C     AT APPROPRIATE RADIUS
      BCKOD=JLOOK
      STNDV=STNDEV(IRAD)
      DEV=BCKOD-ODB
      IF(DEV.LT.0.)DEV=-DEV
      IF(DEV.GT.(3*STNDV))GO TO 9185
C
      JLOOKB=ODB
      IF(JLOOKB.GT.IOVERL)IOVER(J)=1
      IF((JLOOKB+1).LE.IUNDER)IOVER(J)=2
      DJLOOK=ODB-JLOOKB
CHENN>
C      BACKB=LOOKUP(JLOOKB)*(1-DJLOOK)+LOOKUP(JLOOKB+1)*DJLOOK
C      BACK(J)=BACK(J)+LOOKUP(JLOOK)-BACKB
C
      if(JLOOKB.gt.0)then
        BACKB=LOOKUP(JLOOKB)*(1-DJLOOK)+LOOKUP(JLOOKB+1)*DJLOOK
        BACK(J)=BACK(J)+LOOKUP(JLOOK)-BACKB
      else
        BACKB=0.0
      endif
CHENN<
      NGDBCK(J)=NGDBCK(J)+1
      GO TO 185
C     REPLACE WITH VALUE FROM RADIAL DISTRIBUTION CURVE
C     IE BACK(J)=BACK(J)+0
9185  NBDBCK(J)=NBDBCK(J)+1
C
185   CONTINUE
      IF(NFINSH.GT.LINE16) GO TO 200
      JTEST=JBACK(IB)
      JBACK(IB)=MAX0(J+1,JTEST)
      JBDONE(IB)=JBDONE(IB)+1
      JBLEFT(IB)=JBLEFT(IB)-J
      GO TO 200
176   IF(NBEGIN.LE.LINE16) GO TO 174
      IOK=0
      GO TO 200
174   NFIN=NFINSH
      IF(NFINSH.GT.LINE16) GO TO 184
      GO TO 188
184   NFIN = LINE16
188   DO 183 IY=NBEGIN,NFIN
      IL=IY+1-LINE1
      YCO=IY-Y0
      YCOSQ=YCO**2
      DO 183 IX=NXST,NXFI
      JLOOK=JSTRIP(IL,IX)
CHENN>
C      IF(JLOOK.LT.1024.AND.JLOOK.GT.0)GO TO 6996
C      IF(LPRINT) WRITE(6,7002)JH(J),JK(J)
C
      IF(JLOOK.LT.LOOKUPMAX.AND.JLOOK.GT.0)GO TO 6996
      IF(LPRINT) WRITE(6,7002)JH(J),JK(J),JLOOK
CHENN<
      IF(JLOOK.GE.1024)IOVER(J)=1
      IF(JLOOK.LE.1)IOVER(J)=2
      GO TO 183
6996  IF(JLOOK.GT.IOVERL) IOVER(J)=1
      IF(JLOOK.LE.IUNDER) IOVER(J)=2
      XCO=IX-X0
      RAD=1.+SQRT(YCOSQ+XCO**2)
      IRAD=RAD
      IF(RAD.LE.RIN)IOVER(J)=3
      IF (RAD.GT.(ROUT+1.))IOVER(J)=3
      DRAD=RAD-IRAD
      ODB=ODBACK(IRAD)*(1-DRAD)+ODBACK(IRAD+1)*DRAD+YCURVE(IY)
C
      BCKOD=JLOOK
      STNDV=STNDEV(IRAD)
      DEV=BCKOD-ODB
      IF(DEV.LT.0.)DEV=-DEV
      IF(DEV.GT.(3*STNDV))GO TO 9183
C
      JLOOKB=ODB
      IF(JLOOKB.GT.IOVERL)IOVER(J)=1
      IF((JLOOKB+1).LE.IUNDER)IOVER(J)=2
      DJLOOK=ODB-JLOOKB
CHENN>
C      BACKB=LOOKUP(JLOOKB)*(1-DJLOOK)+LOOKUP(JLOOKB+1)*DJLOOK
C      BACK(J)=BACK(J)+LOOKUP(JLOOK)-BACKB
C
      if(JLOOKB.gt.0)then
        BACKB=LOOKUP(JLOOKB)*(1-DJLOOK)+LOOKUP(JLOOKB+1)*DJLOOK
        BACK(J)=BACK(J)+LOOKUP(JLOOK)-BACKB
      else
        BACKB=0.0
      endif
CHENN<
      NGDBCK(J)=NGDBCK(J)+1
      GO TO 183
C
9183  NBDBCK(J)=NBDBCK(J)+1
C
183   CONTINUE
      IF(NFINSH.GT.LINE16) GO TO 200
      JTEST=JBACK(IB)
      JBACK(IB)=MAX0(J+1,JTEST)
      JBDONE(IB)=JBDONE(IB)+1
      JBLEFT(IB)=JBLEFT(IB)-J
200   CONTINUE
      LJSCAN(IB)=JSCAN
      IF(IOK.EQ.0) GO TO 201
      WRITE(6,202)
202   FORMAT(' NOT ENOUGH RANGE IN BACKGROUND SCANS')
      STOP
201   CONTINUE
222   CONTINUE
      LINE1=LINE1+ISLICE
      LINE16=LINE16+ISLICE
      IF(LINE16.GT.NY) GO TO 210
      GO TO 65
210   IF(ISTOP.EQ.0) GO TO 221
      JMEASR=MAX0((JSTART-1),JPASS)
      DO 9209 IB=1,NUMIB
      JINIT=JBACK(IB)
      IF(JINIT.LT.(JBPASS(IB)+1))JINIT=JBPASS(IB)+1
      IF (JINIT.GT.NDATA) GO TO 9209
      DO 209 J=JINIT,NDATA
      IF(LPRINT) WRITE(6,10013)IB,JH(J),JK(J),XBACK(IB,J),YBACK(IB,J)
10013 FORMAT(' FLAG BACKGROUND',I5,' FOR END OF DATA FOR',2I5,
     1' CENTRE',2I10)
      JBDONE(IB)=JBDONE(IB)+1
      JBLEFT(IB)=JBLEFT(IB)-J
209   IFLAG(J)=IFLAG(J)+1
9209  CONTINUE
C
      IF(JDONE.NE.NDATA.OR.JLEFT.NE.0) THEN
        IF(LPRINT) WRITE(6,10018)JLEFT,JDONE
10018 FORMAT(' SOME SPOT MEASUREMENTS MISSING, TOTAL LEFT',I10,
     1' NUMBER OF SPOTS DONE',I10)
      ENDIF
C
      IF(LPRINT) WRITE(6,10017)JMEASR,NDATA,JSKIP
10017 FORMAT(I10,' SPOTS EXAMINED OUT OF',I10,' NUMBER SKIPPED',I5/)
C
C      CHECK THAT ALL BACKGROUND MEASUREMENTS HAVE BEEN DONE
      IF(LPRINT) WRITE(6,10016)
10016 FORMAT(' COMPLETION OF ALL BACKGROUND MEASUREMENTS CHECKED'/)
      DO 9929 IB=1,NUMIB
      IF(JBDONE(IB).EQ.NDATA.AND.JBLEFT(IB).EQ.0)GO TO 9929
      IF(LPRINT) WRITE(6,10015)IB,JBLEFT(IB),JBDONE(IB)
10015 FORMAT(' SOME BACKGROUNDS MISSING FOR',I5,' TOTAL LEFT',I10,
     1' NUMBER OF SPOTS DONE',I10)
9929  CONTINUE
C     DATA HAS NOW BEEN MEASURED
C     REFINE LATTICE PARAMS, REMEASURE, THEN PRINT OUT IF REFINEMENT WAS O.K.
C
      DO 220 J=JBEGIN,JMEASR
      IF(SINT(J).EQ.0.0)GO TO 9230
      PROP=IPROP(J)/1000.
      NXMRA=PROP*(NXMT-NXM)+NXM
      NYMRA=PROP*(NYMT-NYM)+NYM
      NXBRA=PROP*(NXBT-NXB)+NXB
      NYBRA=PROP*(NYBT-NYB)+NYB
      AMEAS=NXMRA*NYMRA
      ABACK=NUMIB*NXBRA*NYBRA
      FACTOR=AMEAS/ABACK
      IF (IFLAG(J).NE.0)GO TO 225
      ICORR(J)=ICORR(J)-BACK(J)*FACTOR
      GO TO 220
C    INSIST ON AT LEAST 3 BACKGROUND RASTERS OUT OF 6; OR 2 OUT OF 4 
C    IF NTYPE=2
225   IF(IFLAG(J).GT.2.AND.NTYPE.EQ.2)GO TO 230 
      IF(IFLAG(J).GT.4) GO TO 230
      ABACK=(NUMIB-IFLAG(J))*NXBRA*NYBRA
      FSPESL=AMEAS/ABACK
      ICORR(J)=ICORR(J)-BACK(J)*FSPESL
      IF(LPRINT) WRITE(6,226) JH(J),JK(J),ICORR(J),IFLAG(J)
226   FORMAT(2I5,I10,'  THIS SPOT CORR. FOR BACKGROUND WITH',I5,
     1'  MISSING BACKGROUND RASTERS ')
      GO TO 220
9230  IF(LPRINT) WRITE(6,9232)JH(J),JK(J)
9232  FORMAT(2I5,' THIS SPOT DELETED BECAUSE ZERO SPOT RASTER')
      ICORR(J)=0
      GO TO 220
230   ICORR(J)=0
      SINT(J)=0.0
      IF(LPRINT) WRITE(6,231) JH(J),JK(J),IFLAG(J)
231   FORMAT(2I5,'  THIS SPOT DELETED BECAUSE OF',I5,'  MISSING BACKGROU
     1ND RASTERS OUT OF SIX')
220   CONTINUE
      IF(ISTOP.EQ.1) GO TO 300
221   CALL SEARCH(X0NEW,Y0NEW,DX1NEW,DY1NEW,DX2NEW,DY2NEW,JH,JK,
     1XCOREF,YCOREF,SINT,NDATA,B3,XEWALD,YEWALD,X0,Y0,
     2DX1,DY1,DX2,DY2)
C
      WRITE(6,9211)
9211  FORMAT(' REFINED VALUES FOR INPUT TO FURTHER CYCLES ARE THOSE ',
     1'IN LINE ABOVE'/)
      WRITE(6,211) X0NEW,Y0NEW
      WRITE(6,212) DX1NEW,DY1NEW,DX2NEW,DY2NEW
211   FORMAT('  REFINED VALUES OF X0 AND Y0  ',2F10.2)
212   FORMAT('  REFINED VECTOR FOR  (1,0)  ',2F10.3/'  REFINED VECTOR
     1FOR (0,1)  ',2F10.3)
      IF(ABS(X0NEW-X0).GT.0.2) GO TO 250
      IF(ABS(Y0NEW-Y0).GT.0.2) GO TO 250
      IF(ABS(DX1NEW-DX1).GT.0.05) GO TO 250
      IF(ABS(DY1NEW-DY1).GT.0.05) GO TO 250
      IF(ABS(DX2NEW-DX2).GT.0.05) GO TO 250
      IF(ABS(DY2NEW-DY2).GT.0.05) GO TO 250
C     AT LEAST 2 CYCLES TO BE DONE IF NOT PRESET
      IF(NCYCLE.LT.2.AND.NCYC.NE.1) GO TO 250
      WRITE(6,213) NCYCLE
213   FORMAT(' END OF LATTICE PARAM REFINEMENT  '/' NUMBERS CONVERGE',
     1' AFTER',I5,' CYCLES')
9252      MDONE=1
      AH=DX1NEW**2+DY1NEW**2
      AK=DX2NEW**2+DY2NEW**2
      AI=(DX2NEW-DX1NEW)**2+(DY2NEW-DY1NEW)**2
      AH=SQRT(AH)
      AK=SQRT(AK)
      AI=SQRT(AI)
C
C  THESE ANGLES DEPEND ON SCANNING AXIS CONVENTION WHICH AT JUNE 82
C  IS ORIGIN IN TOP RIGHT CORNER NEAR THE FILM NUMBER.
      ANGH=ATAN2(DY1NEW,DX1NEW)
      ANGK=ATAN2(DY2NEW,DX2NEW)
      ANGI=ATAN2(-DY1NEW+DY2NEW,DX2NEW-DX1NEW)
      ANGH=ANGH*180.0/PI
      ANGK=ANGK*180.0/PI
      ANGI=ANGI*180.0/PI
      WRITE(6,214)AH,AK,AI
      WRITE(6,215)ANGH,ANGK,ANGI
214   FORMAT(60X,'LENGTHS OF H,K,-H+K  ARE',3F10.3)
215   FORMAT(60X,'ANGLES TO PLATE EDGE ARE',3(F7.0,3X))
C
C  CALCULATION OF EXACT TLTAXA,TLTANG FROM LATTICE PARAMETERS & DIST CORR.
      COSDIS=COS(ANGDIS)
      SINDIS=SIN(ANGDIS)
      DX1COR=COSDIS*(DX1NEW*COSDIS+DY1NEW*SINDIS)-SHRINK*SINDIS*
     1(-DX1NEW*SINDIS+DY1NEW*COSDIS)
      DY1COR=SINDIS*(DX1NEW*COSDIS+DY1NEW*SINDIS)+SHRINK*COSDIS*
     1(-DX1NEW*SINDIS+DY1NEW*COSDIS)
      DX2COR=COSDIS*(DX2NEW*COSDIS+DY2NEW*SINDIS)-SHRINK*SINDIS*
     1(-DX2NEW*SINDIS+DY2NEW*COSDIS)
      DY2COR=SINDIS*(DX2NEW*COSDIS+DY2NEW*SINDIS)+SHRINK*COSDIS*
     1(-DX2NEW*SINDIS+DY2NEW*COSDIS)
      WRITE(6,9987)DX1COR,DY1COR,DX2COR,DY2COR
9987  FORMAT(' DISTORTION CORRECTED REFINED VECTORS, (1,0)',2F10.3/
     1       '                                       (0,1)',2F10.3)
      PI=3.141592654
      ASTILT=SQRT(DX1COR**2+DY1COR**2)
      BSTILT=SQRT(DX2COR**2+DY2COR**2)
      ANGA=ATAN2(DY1COR,DX1COR)
      ANGB=ATAN2(DY2COR,DX2COR)
      DANG=ANGA-ANGB
      IF(DANG.GT.PI)DANG=DANG-2.0*PI
      IF(DANG.LT.-PI)DANG=DANG+2.0*PI
      GSTILT=DANG*180.0/PI
      IF(GSTILT.LT.0.0) GSTILT=-GSTILT
      CALL EMTILT(TL,TAXA,TANG,ASTAR,BSTAR,GMSTAR,ASTILT,BSTILT,GSTILT)
      WRITE(6,8004)TAXA,TANG,TL
C
      ISTOP =1
      GO TO 9251
9250  WRITE(6,9213)NCYCLE
9213  FORMAT(' END OF FIXED NUMBER OF CYCLES',I5)
      ISTOP=1
      GO TO 9252
250   IF(NCYCLE.EQ.NCYC.AND.MDONE.NE.1)GO TO 9250
9251  X0=X0NEW
      Y0=Y0NEW
      DX1=DX1NEW
      DY1=DY1NEW
      DX2=DX2NEW
      DY2=DY2NEW
      WRITE(6,251) NCYCLE
251   FORMAT(' END OF CYCLE NUMBER  ',I5,'   ***************************
     1***********************************************************')
      NCYCLE=NCYCLE +1
C
      CALL IMPOSN(1,0,0)
      GO TO 260
300   CONTINUE
C     NOW OUTPUT OF ALL DATA
C     FIRST SIMPLY A LIST OF ALL NUMBERS  TOGETHER WITH CENTRE DEVIATIONS
      IF(LPRINT) WRITE(6,399)
      IF(LPRINT) WRITE(6,312)
      IF(LPRINT) WRITE(6,9312)
C     STATISTICS OF CENTRE OF GRAVITY DISPLACEMENTS
C     CALCULATED AS A FUNCTION OF RADIUS IN PATTERN FROM XD,YD
C
C     PLOT C OF G DEVIATIONS
      PLTSIZ = 260.0      ! old plots were dimensioned in mm
      FONTSIZE = 4.      ! set fontsize to 4mm
      CALL P2K_OUTFILE('PLOTOUT.PS',10)
      CALL P2K_HOME
      CALL P2K_FONT('Courier'//CHAR(0),FONTSIZE)
      CALL P2K_GRID(0.5*PLTSIZ,0.5*PLTSIZ,1.0)
      CALL P2K_ORIGIN(-0.5*PLTSIZ,-0.7*PLTSIZ,0.)
      CALL P2K_COLOUR(0)
      CALL P2K_MOVE(0.,0.,0.)
      XWIDTH = PLTSIZ*FLOAT(NX)/MAX0(NX,NY)
      YHEIGHT = PLTSIZ*FLOAT(NY)/MAX0(NX,NY)
      CALL P2K_DRAW(XWIDTH,0.,0.)
      CALL P2K_DRAW(XWIDTH,YHEIGHT,0.)
      CALL P2K_DRAW(0.,YHEIGHT,0.)
      CALL P2K_DRAW(0.,0.,0.)
C
      CALL P2K_MOVE(XWIDTH-20.0,YHEIGHT,0.)
      CALL P2K_DRAW(XWIDTH-20.0,YHEIGHT-10.0,0.)
      CALL P2K_DRAW(XWIDTH,YHEIGHT-10.0,0.)
      CALL P2K_MOVE(XWIDTH-21.5,YHEIGHT-6.5,0.)
CTSH      WRITE(TEXT,9300)NPLATE
CTSH++
      WRITE(TEXT,9300)NPLATE
CTSH--
9300  FORMAT(I10)
      CALL P2K_STRING(TEXT,6,0.)
      CALL P2K_MOVE(XWIDTH/2.0,0.,0.)
      CALL P2K_COLOUR(0)
      CALL P2K_DRAW(XWIDTH/2.0,YHEIGHT,0.)
      CALL P2K_MOVE(0.,YHEIGHT/2.0,0.)
      CALL P2K_DRAW(XWIDTH,YHEIGHT/2.,0.)
C      SPLOT=PLTSIZ/(MAX0(NX,NY)-1.0)
      SPLOT=PLTSIZ/MAX0(NX,NY)
      XOFFSET = 0.
      YOFFSET = -0.8
      X0PLOT=XWIDTH-X0*SPLOT + XOFFSET
      Y0PLOT=YHEIGHT-Y0*SPLOT + YOFFSET
      CALL P2K_MOVE(X0PLOT,Y0PLOT,0.)
CTSH      ENCODE (1,9301,TEXT)
CTSH++
      WRITE (TEXT,9301)
CTSH--
      CALL P2K_CSTRING(TEXT,1,0.)
9301  FORMAT('*')
      CALL P2K_COLOUR(0)
      CALL P2K_FONT('Courier'//CHAR(0),0.5*FONTSIZE)

C
      RDEL = ( ROUT )/8.0
      RMST=0.0
      RT=0.0
      NT=0
      INTMAX=0
      DO 307 J=1,8
      DELRAD(J)=0.0
      THETA(J)=0.0
307   NRAD(J)=0
      DO 308 J=JBEGIN,JMEASR
      PROP=IPROP(J)/1000.
      NXBRA=PROP*(NXBT-NXB)+NXB
      NYBRA=PROP*(NYBT-NYB)+NYB
      ABACK=NUMIB*NXBRA*NYBRA
      BACK(J)=BACK(J)/ABACK
      NXMRA=PROP*(NXMT-NXM)+NXM
      NYMRA=PROP*(NYMT-NYM)+NYM
C
      IF(INTMAX.LT.ICORR(J).AND.IOVER(J).EQ.0) INTMAX=ICORR(J)
      IF(XCOREF(J).EQ.0.0) GO TO 309
      CGX=XCOREF(J)-XCOORD(J)
      CGY=YCOREF(J)-YCOORD(J)
      RMST=RMST+(CGX*CGX+CGY*CGY)
      RT=RT+SQRT(CGX*CGX+CGY*CGY)
      NT=NT+1
      XVECT=XCOORD(J)-X0
      YVECT=YCOORD(J)-Y0
      RVECT=SQRT(XVECT*XVECT+YVECT*YVECT)
      IRSLOT=RVECT/RDEL+1
      IF(IRSLOT.GT.8) IRSLOT=8
      COSTH=XVECT/RVECT
      SINTH=YVECT/RVECT
      DELRAD(IRSLOT)=DELRAD(IRSLOT)+CGX*COSTH+CGY*SINTH
      THETA(IRSLOT)=THETA(IRSLOT)+CGX*SINTH-CGY*COSTH
      NRAD(IRSLOT)=NRAD(IRSLOT)+1
C
      GO TO 310
309   CONTINUE
      XPLOT=XWIDTH-XCOORD(J)*SPLOT + XOFFSET
      YPLOT=YHEIGHT-YCOORD(J)*SPLOT + YOFFSET
      CALL P2K_MOVE(XPLOT,YPLOT,0.)
CTSH      IF(IOVER(J).EQ.0) WRITE(TEXT,9303)
CTSH      IF(IOVER(J).EQ.2) WRITE(TEXT,9304)
CTSH      IF(IOVER(J).EQ.1) WRITE(TEXT,9305)
CTSH      IF(IOVER(J).EQ.3) WRITE(TEXT,9306)
CTSH++
      IF(IOVER(J).EQ.0) WRITE(TEXT,9303)
      IF(IOVER(J).EQ.2) WRITE(TEXT,9304)
      IF(IOVER(J).EQ.1) WRITE(TEXT,9305)
      IF(IOVER(J).EQ.3) WRITE(TEXT,9306)
CTSH--
      CALL P2K_CSTRING(TEXT,1,0.)
9303  FORMAT('O')
9304  FORMAT('-')
9305  FORMAT('+')
9306  FORMAT('R')
      IF(IOVER(J).NE.0)GO TO 9309
      IF(SINT(J).EQ.0.0)GO TO 9309
      IF(IFLAG(J).GT.2)GO TO 9309
C
      IF(LPRINT) WRITE(6,311) JH(J),JK(J),ICORR(J),IPROP(J),
     $NXMRA,NYMRA,SINT(J),BACK(J),NXBRA,NYBRA,NBDBCK(J),NGDBCK(J)
      GO TO 308
C  
9309  IF(LPRINT) WRITE(6,9311)JH(J),JK(J),ICORR(J),IPROP(J),
     $NXMRA,NYMRA,SINT(J),BACK(J),NXBRA,NYBRA,NBDBCK(J),NGDBCK(J)
9311  FORMAT(2I5,I8,I7,1X,'(',I2,',',I2,')**',F10.0,F12.1,
     15X,'(',I2,',',I2,')',I5,I6)
      GO TO 308
310   IF(LPRINT) WRITE(6,311) JH(J),JK(J),ICORR(J),IPROP(J),NXMRA,
     $NYMRA,SINT(J),BACK(J),NXBRA,NYBRA,NBDBCK(J),NGDBCK(J),
     $CGX,CGY,XVECT,YVECT
      XPLOT=XWIDTH-XCOORD(J)*SPLOT + XOFFSET
      YPLOT=YHEIGHT-YCOORD(J)*SPLOT + YOFFSET
      CALL P2K_MOVE(XPLOT,YPLOT,0.)
CTSH      WRITE(TEXT,9302)
CTSH++
      WRITE(TEXT,9302)
CTSH--
9302  FORMAT('X')
      CALL P2K_CSTRING(TEXT,1,0.)
      XPLOT=XWIDTH-XCOORD(J)*SPLOT
      YPLOT=YHEIGHT-YCOORD(J)*SPLOT
      CALL P2K_MOVE(XPLOT,YPLOT,0.)
      XPLOT=XWIDTH-SPLOT*(20*CGX+XCOORD(J))
      YPLOT=YHEIGHT-SPLOT*(20*CGY+YCOORD(J))
      IF(XPLOT.GT.XWIDTH.OR.XPLOT.LT.0.)GO TO 308
      IF(YPLOT.GT.YHEIGHT.OR.YPLOT.LT.0.)GO TO 308
      CALL P2K_DRAW(XPLOT,YPLOT,0.)
308   CONTINUE
      CALL P2K_PAGE
311   FORMAT(2I5,I8,I7,1X,'(',I2,',',I2,')',2X,F10.0,F12.1,
     15X,'(',I2,',',I2,')',I5,I6,9X,2F6.2,5X,2F10.2)
312   FORMAT('   IH   IK     INT   PROP RAST DIM      PEAK      AVRGE (B
     1CK  BCK RAST NBAD NGOOD         CGX   CGY          XVECT     Y',
     2'VECT'/)
9312  FORMAT('             (CORR)                  (TOTAL OD)   -RADL CO
     1MP)   DIM'//)
C     OUTPUT STATISTICS OF CENTRE OF GRAVITY
      WRITE(6,321)
      DO 320 J=1,8
      IF(NRAD(J).EQ.0) GO TO 318
      THETA(J)=THETA(J)/NRAD(J)
      DELRAD(J)=DELRAD(J)/NRAD(J)
318   ROUT=RDEL*J
      RCORR = B3*ROUT*ROUT*ROUT
      WRITE(6,322) ROUT,DELRAD(J),THETA(J),NRAD(J),RCORR
320   CONTINUE
      IF(NT.EQ.0)GO TO 6324
      RT=RT/NT
      RMST=SQRT(RMST/NT)
6324  WRITE(6,323) RT,NT
      WRITE(6,324) RMST,NT
324   FORMAT(' OVERALL RMS DEVIATION=',F6.3,'  FOR TOTAL NUMBER=',I5)
323   FORMAT(' OVERALL AVERAGE DEVTN=',F6.3,'  FOR TOTAL NUMBER=',I5)
321   FORMAT(/' DEVIATION OF CENTRE OF GRAVITY FROM PREDICTED POSITION'/
     1'  RADIUS   RADIAL  ANGULAR  NUMBER OF SPOTS   APPROX RADIAL'/
     2'             COMPONENTS                       CORR APPLIED'/)
322   FORMAT(F8.1,2F9.3,I8,12X,F9.3)
C
C
      CALL STATS(NDATA,JH,JK,ICORR,IOVER,
     .		X0,Y0,XCOORD,YCOORD,ASTAR,AHCORR,ROUTSTATS,RFA,NTOT)

C
C     NOW AVERAGE SYMMETRY RELATED SPOTS(ONLY FRIEDEL RELATION).
C     AND OUTPUT AVERAGED DATA
      IF(LPRINT) WRITE(6,399)
      NREJECTOUR=0
      NREJECTO=0
      NREJECTU=0
      NREJECTR=0
      NMISS=0
      NDISAGREE=0
      TOP=0.0
      BOT=0.0
      NNEG=0
      NPOS=0
      JSPOTS=0
      JAVSTP=INTMAX/10
      JDISTP=INTMAX/100
      DO 350 I1=1,10
      TOPT(I1)=0.0
      BOTT(I1)=0.0
      DO 350 I2=1,12
350   ITABLE(I1,I2)=0
      IF(LPRINT) WRITE(6,351)
351   FORMAT(' FRIEDEL PAIRS NOW AVERAGED'//
     1'   IH   IK       JAV     JDIFF           AV/PNT',
     2'   DIF/PNT       RASTER AREAS')
      IF(LPRINT) WRITE(6,355) FRACT,ABSOL
355   FORMAT(45X,' REJECTED IF JDIFF IS GREATER THAN',F5.2,' *JAV, AND',
     1F10.2)
      WRITE(6,404) NPLATE,TITLE
      WRITE(2,406) NPLATE,TITLE
      DO 9402 J=JBEGIN,JMEASR
      IF(JH(J).EQ.0.AND.JK(J).EQ.0)GO TO 9402
      IMATCH=0
      DO 400 K=J,JMEASR
      IF(JH(J).NE.-JH(K)) GO TO 400
      IF(JK(J).NE.-JK(K)) GO TO 400
      IF((SINT(J).EQ.0.).OR.(SINT(K).EQ.0.)) GO TO 360
      JDIFF = ICORR(J)-ICORR(K)
      JAV=(ICORR(J)+ICORR(K))/2
      MDIFF=JDIFF
      IF(JDIFF.LT.0) MDIFF=-JDIFF
      JFRACT=FRACT*JAV
      JABSOL=ABSOL
      PROP=IPROP(J)/1000.
      NXMRAJ=PROP*(NXMT-NXM)+NXM
      NYMRAJ=PROP*(NYMT-NYM)+NYM
      AMEASJ=NXMRAJ*NYMRAJ
      PROP=IPROP(K)/1000.
      NXMRAK=PROP*(NXMT-NXM)+NXM
      NYMRAK=PROP*(NYMT-NYM)+NYM
      AMEASK=NXMRAK*NYMRAK
C
      AVPNT=0.5*(ICORR(J)/AMEASJ+ICORR(K)/AMEASK)
      DFPNT=ICORR(J)/AMEASJ-ICORR(K)/AMEASK
      IF((IOVER(J).EQ.1).OR.(IOVER(K).EQ.1)) GO TO 356
      IF((IOVER(J).EQ.2).OR.(IOVER(K).EQ.2)) GO TO 356
      IF((IOVER(J).EQ.3).OR.(IOVER(K).EQ.3)) GO TO 356
      IF((MDIFF.GT.JFRACT).AND.(MDIFF.GT.JABSOL)) GO TO 352
      GO TO 353
360   IF(LPRINT) WRITE(6,361)JH(J),JK(J),ICORR(J),ICORR(K)
361   FORMAT(2I5,20X,'SPOT MISSING, INTENSITIES WERE  ',2I10)
      NMISS=NMISS+1
      GO TO 9404
356   IF(LPRINT) WRITE(6,357) JH(J),JK(J),JAV,JDIFF,IOVER(J),IOVER(K)
357   FORMAT(2I5,2I10,'     REFLECTION REJECTED, OPTICAL DENSITY TOO ',
     1' HIGH OR LOW ON SPOT OR BACKGROUND',2I5)
      	IF((IOVER(J).NE.0).OR.(IOVER(K).NE.0))NREJECTOUR=NREJECTOUR+1
      	IF((IOVER(J).EQ.1).OR.(IOVER(K).EQ.1))NREJECTO=NREJECTO+1
      	IF((IOVER(J).EQ.2).OR.(IOVER(K).EQ.2))NREJECTU=NREJECTU+1
      	IF((IOVER(J).EQ.3).OR.(IOVER(K).EQ.3))NREJECTR=NREJECTR+1
      GO TO 9404
C
352   IF(LPRINT) WRITE(6,354) JH(J),JK(J),JAV,JDIFF,AVPNT,DFPNT,
     $NXMRAJ,NYMRAJ,NXMRAK,NYMRAK
354   FORMAT(2I5,2I10,5X,2F10.1,2(7X,I2,',',I2),'   REFLECTION REJ',
     1'ECTED  JDIFF TOO BIG')
      NDISAGREE=NDISAGREE+1
      GO TO 9404
353   JSPOTS=JSPOTS+1
      IF(LPRINT) WRITE(6,401) JH(J),JK(J),JAV,JDIFF,AVPNT,DFPNT,
     $NXMRAJ,NYMRAJ,NXMRAK,NYMRAK
      IF(JAV.LE.0) NNEG=NNEG+1
      IF(JAV.GT.0) NPOS=NPOS+1
      IF(JAV.EQ.0) JAV=1
      IF(JAV.LT.MININT) JAV=MININT
      JSAV=JAV
      IF(JSAV.LE.0)JSAV=1
      TOP=TOP+MDIFF
      BOT=BOT+JSAV
      IAVSLT=(JSAV/JAVSTP)+2
      IF(IAVSLT.EQ.2) IAVSLT=((JSAV*2)/JAVSTP)+1
      IDISLT=(MDIFF/JDISTP)+1
      IF(IAVSLT.GT.10) IAVSLT=10
      IF(IDISLT.GT.12) IDISLT=12
      IF(IAVSLT.LT.1) IAVSLT=1
      IF(IDISLT.LT.1)IDISLT=1
      ITABLE(IAVSLT,IDISLT)=ITABLE(IAVSLT,IDISLT)+1
      TOPT(IAVSLT)=TOPT(IAVSLT)+MDIFF
      BOTT(IAVSLT)=BOTT(IAVSLT)+JSAV
CHENN>
CAND MDIFF instead of JDIFF because JDIFF is sometimes negative.
C      WRITE(2,405) JH(J),JK(J),JAV,JDIFF
      WRITE(2,405) JH(J),JK(J),JAV,MDIFF
CHENN<
C
C SET INDICES TO ZERO AFTER REFLECTIONS DEALT WITH
9404  JH(K)=0
      JK(K)=0
      IMATCH=1
      GO TO 9403
400   CONTINUE
9403  IF(IMATCH.EQ.0)GO TO 9402
      JH(J)=0
      JK(J)=0
9402  CONTINUE
      JHTERM=100
      WRITE(2,405) JHTERM,JHTERM,JHTERM,JHTERM
      IF(LPRINT) WRITE(6,401) JHTERM,JHTERM,JHTERM
C
C    SEARCH FOR ANY UNPAIRED REFLECTIONS
      NOPAIR=0
      WRITE(6,9351)
9351  FORMAT(//' REFLECTIONS WITHOUT FRIEDEL PAIRS'//
     1'   IH   IK     ICORR'//)
      DO 9400 J=JBEGIN,JMEASR
      IF (JH(J).EQ.0.AND.JK(J).EQ.0)GO TO 9400
      NOPAIR=NOPAIR+1
      WRITE(6,9401)JH(J),JK(J),ICORR(J)
9401  FORMAT(2I5,I10)
9400  CONTINUE
C
C     OUTPUT OF STATISTICAL TABLES AND R-FACTORS
C
      WRITE(6,506)
      WRITE(6,500)JDISTP
      DO 520 K=1,12
520   ITAV(K)=K*JDISTP
      WRITE(6,501) (ITAV(K),K=1,12)
      WRITE(6,502)
      ITOTAL=0
      DO 550 I1=1,10
      IALL=0
      IF(I1.GT.2) GO TO 522
      IAVB=((I1-1)*JAVSTP)/2
      IAVF=(I1*JAVSTP)/2
      GO TO 523
522   IAVB=(I1-2)*JAVSTP
      IAVF=(I1-1)*JAVSTP
523   CONTINUE
      DO 540 I2=1,12
540   IALL=IALL+ITABLE(I1,I2)
      IF(IALL.EQ.0) GO TO 550
      SYM=TOPT(I1)/BOTT(I1)
      WRITE(6,503)IAVB,IAVF,(ITABLE(I1,K),K=1,12),IALL,SYM
      ITOTAL=ITOTAL+IALL
550   CONTINUE
      WRITE(6,504)
      DO 560 I2=1,12
      DO 560 I1=2,10
560   ITABLE(1,I2)=ITABLE(1,I2) + ITABLE(I1,I2)
      WRITE(6,505)(ITABLE(1,K),K=1,12),ITOTAL
      WRITE(6,507) NPOS,NNEG
      WRITE(6,399)
      RSYM=TOP/BOT
      WRITE(6,402) RSYM,JSPOTS
      IF(NREJECTOUR.NE.0) WRITE(6,4407) NREJECTOUR
      IF(NREJECTO.NE.0) WRITE(6,4401) NREJECTO
      IF(NREJECTU.NE.0) WRITE(6,4406) NREJECTU
      IF(NREJECTR.NE.0) WRITE(6,4402) NREJECTR
      IF(NDISAGREE.NE.0) WRITE(6,4403) NDISAGREE
      IF(NOPAIR.NE.0) WRITE(6,4404) NOPAIR
      IF(NMISS.NE.0) WRITE(6,4405) NMISS
      WRITE(6,4400) RFA,NTOT
      STOP
513   FORMAT(/' PICKAUTO - VX4.00 (18.8.00) - automatic ',
     .	'electron diffraction spot integration'/) 
512   FORMAT('  RADIAL DENSITY OR Y-CURVE FOR WRONG FILM ON DISC')
511   FORMAT(I10,2F11.5)
510   FORMAT(//' RADIAL DENSITY'/I10,18A4)
509   FORMAT(2F11.5)
508   FORMAT(3F10.2,2I10)
507   FORMAT(//' TOTAL POSITIVE SPOTS=',I5/' TOTAL NEGATIVE SPOTS=',I5)
506   FORMAT(///' OUTPUT OF STATISTICS AND R-FACTORS'///)
505   FORMAT(/24X,'ALL  ',12I5,I10)
504   FORMAT(18X,'AND ABOVE')
500   FORMAT(29X,'JDIFF IN GROUPS OF',I6/)
501   FORMAT(24X,'JAV  ',12I5,'       ALL    RSYM')
502   FORMAT(82X,'AND ABOVE'/)
503   FORMAT(5X,I5,5X,'TO',5X,I5,2X,12I5,I10,F9.3)
406   FORMAT(I10,18A4,'***')
405   FORMAT(2I5,2I6)
404   FORMAT(5X,I10,18A4)
CHENN>
C403   FORMAT(I5,18A4)
403   FORMAT(I10,18A4)
CHENN<
402   FORMAT('  OVERALL MEAN DIFFERENCE BET SYM REL SPOTS =  ',F7.3/
     .'  FOR  ',I5,'  SPOTS, with the following not written out')
4400   FORMAT(/'  OVERALL MEAN DIFFERENCE BET SYM REL SPOTS =  ',F7.3/
     .'  FOR  ',I5,'  SPOTS, with only overloads excluded')
4401  FORMAT(25X,I5,'  pairs rejected for overloads')
4402  FORMAT(25X,I5,'  pairs rejected beyond radius limits')
4403  FORMAT(20X,I5,'  pairs rejected for disagreement (FRACT,ABSOL)')
4404  FORMAT(20X,I5,'  individual spots with no Friedel pair')
4405  FORMAT(20X,I5,'  unexplained missing pairs')
4406  FORMAT(25X,I5,'  pairs rejected for underloads(pointer?)')
4407  FORMAT(20X,I5,'  pairs rejected for all reasons (O,U,R)')
401   FORMAT(2I5,2I10,5X,2F10.1,2(7X,I2,',',I2),'   is OK')
399   FORMAT(/////)
777   FORMAT(' SMALL RASTER BEFORE BIG RASTER IN LIST, FUNNY!')
C
8999  WRITE(6,8998)
8998  FORMAT(' ERROR ON INPUT FILE')
      STOP
C

590   WRITE(6,512)
      STOP
      END
C
C##############################################################################
C
      SUBROUTINE CENTRE(JSTORE,IS,XCEN,YCEN,ABNORM,IPRINT,NXM,NYM,JH,JK,
     1XAMINE,ODBACK,ODBTOT,NXST,NYST,X0,Y0,
     2ROUT,RIN,YCURVE,
     3LOOKUP,IOVERL,IUNDER,IOVERC,ISTOP)
C
CHENN>
C      PARAMETER (NMAX=6000)
C
      PARAMETER (NMAX=48000)
      PARAMETER (LOOKUPMAX=10240)
CHENN<
      PARAMETER (ISPOTS=140)
      PARAMETER (IRAST=80)
      INTEGER*2 JSTORE(ISPOTS,IRAST,IRAST),JH,JK
      DIMENSION ODBACK(1),YCURVE(1)
      REAL*4 LOOKUP(LOOKUPMAX)
C  check for overloads and underloads on all numbers 12-Nov-1994
C
C     THIS SUBROUTINE CALCULATES THE CENTRE OF GRAVITY OF EACH SPOT
C     IT TERMINATES ABNORMALLY (ABNORM=1,0) IF THE SPOT IS TOO WEAK
C     OR TOO STRONG, AND PRINTS OUT SOME DATA IF IPRINT=1
C     NEW VERSION 19.5.80 INCORPORATES RADIAL DENSITY SUBTRACTION.
C
C     CRITERION FOR DETECTION OF A SPOT IS THAT THERE SHOULD BE
C     MORE THAN 10 PERCENT OF THE MEASURED RASTER AREA HAVING A
C     DENSITY GREATER THAN IMIN + XAMINE*SQRT(MEAN PERIMETER OD).
C     IMIN IS LOWEST AVERAGE DENSITY ALONG TWO ADJACENT EDGES OF RASTER.
C
C    AFTER FINAL CYCLE ISTOP=1; CALCULATE TOTAL RADIAL AND Y-CURVE BACKGROUND
C    CORRECTION AND RETURN TO MAIN AS ODBTOT
C    SET OVERLOAD INDICATORS IF NECESSARY
C
C    DO NOT BOTHER WITH Y-CURVE CORRECTION IN C.OF G. AND SPOT DETECTION
C    PARTS, BUT DO INCLUDE IN FINAL CYCLE.
C
      X=XAMINE
      IF(X.EQ.0.0)X=1.5
      IF(ABNORM.LT.0.0)GO TO 1000
      ABNORM=0.0
1000  IMAX=0
      IMIN=IOVERL
      TOT=0.
      PER1=0.
      PER2=0.
      PER3=0.
      PER4=0.
      FACTOR=100./(NXM*NYM)
      LIMIT=(NXM*NYM)/10
      X1=NXST-X0
      X2=X1+NXM-1
      Y1=NYST-Y0
      Y2=Y1+NYM-1
      X1SQ=X1**2
      X2SQ=X2**2
      Y1SQ=Y1**2
      Y2SQ=Y2**2
      DO 5 IX=1,NXM
      TOT=TOT+JSTORE(IS,IX,1)+JSTORE(IS,IX,NYM)
      XCOORD=IX-1+NXST-X0
      XCOSQ=XCOORD**2
      IRAD=1.5+SQRT(XCOSQ+Y1SQ)
      IF(IRAD.GT.(ROUT+1))GO TO 950
      IF(IRAD.LE.RIN)GO TO 950
      PER1=PER1+JSTORE(IS,IX,1)-ODBACK(IRAD)
      IRAD=1.5+SQRT(XCOSQ+Y2SQ)
      IF(IRAD.GT.(ROUT+1))GO TO 950
      IF(IRAD.LE.RIN)GO TO 950
5     PER2=PER2+JSTORE(IS,IX,NYM)-ODBACK(IRAD)
      PER1=PER1/NXM
      PER2=PER2/NXM
      DO 10 IY=1,NYM
      TOT=TOT+JSTORE(IS,1,IY)+JSTORE(IS,NXM,IY)
      YCOORD=IY-1+NYST-Y0
      YCOSQ=YCOORD**2
      IRAD=1.5+SQRT(YCOSQ+X1SQ)
      IF(IRAD.GT.(ROUT+1))GO TO 950
      IF(IRAD.LE.RIN)GO TO 950
      PER3=PER3+JSTORE(IS,1,IY)-ODBACK(IRAD)
      IRAD=1.5+SQRT(YCOSQ+X2SQ)
      IF(IRAD.GT.(ROUT+1))GO TO 950
      IF(IRAD.LE.RIN)GO TO 950
10    PER4=PER4+JSTORE(IS,NXM,IY)-ODBACK(IRAD)
      PER3=PER3/NYM
      PER4=PER4/NYM
      TOT=TOT/(2*(NXM+NYM))
      AV1=(PER1+PER3)/2.
      AV2=(PER2+PER3)/2.
      AV3=(PER2+PER4)/2.
      AV4=(PER4+PER1)/2.
      XMIN=AMIN1(AV1,AV2,AV3,AV4)
      IMIN=XMIN
C     NOW CALCULATE CENTRE OF GRAVITY
      ITEST=XMIN+X*SQRT(TOT)
      NABOVE=0
      DO 15 IX=1,NXM
      XCOORD=NXST+IX-1-X0
      XCOSQ=XCOORD**2
      DO 15 IY=1,NYM
      YCOORD=NYST+IY-1-Y0
      IRAD=1.5+SQRT(XCOSQ+YCOORD**2)
      IF(JSTORE(IS,IX,IY).GT.IMAX) IMAX=JSTORE(IS,IX,IY)
      IF(JSTORE(IS,IX,IY).LT.IMIN) IMIN=JSTORE(IS,IX,IY)
      IF((JSTORE(IS,IX,IY)-ODBACK(IRAD)).GE.ITEST) NABOVE=NABOVE+1
15    CONTINUE
      IF(NABOVE.LE.LIMIT) GO TO 50
CHENN>
C      IF(IMIN.LT.IUNDER.OR.IMAX.GT.750) THEN
C
      IF(IMIN.LT.IUNDER.OR.IMAX.GT.IOVERL) THEN
CHENN<
      	IF(IPRINT.NE.0) THEN
CHENN>
C      	  WRITE(6,16) JH,JK,IMIN,IMAX
C16	  FORMAT(' spot',2I5,' overload or underload, IMIN<1,IMAX>750',2I5)
C
          WRITE(6,16) JH,JK,IUNDER,IOVERL,IMIN,IMAX
16        FORMAT(' spot',2I5,' overload or underload, IMIN<',I5,
     .     ',IMAX>',I5,2I5)
CHENN<
      	  IXSTP=1
      	  IYSTP=1
      	  IF(NXM.LE.30)GO TO 2012
      	  IXSTP=2
      	  IYSTP=2
2012  	  DO 2014 IY=1,NYM,IYSTP
2014      WRITE(6,12)(JSTORE(IS,IX,IY),IX=1,NXM,IXSTP)
      	ENDIF
      	GO TO 50
      ENDIF
C		  overload test different from IOVERL
      ISUM=0
      IXSUM=0
      IYSUM=0
      DO 30 IX=1,NXM
      XCOORD=NXST+IX-1-X0
      XCOSQ=XCOORD**2
      DO 30 IY=1,NYM
      YCOORD=NYST+IY-1-Y0
      IRAD=1.5+SQRT(XCOSQ+YCOORD**2)
      ID=JSTORE(IS,IX,IY)-ODBACK(IRAD)-XMIN
      ISUM=ISUM+ID
      IXSUM=IXSUM+ID*IX
30    IYSUM=IYSUM+ID*IY
      XSUM=IXSUM
      YSUM=IYSUM
      SUM=ISUM
      IF(SUM.NE.0.0)GO TO 930
      XCEN=NXM/2.
      YCEN=NYM/2.
      GO TO 931
930   XCEN=XSUM/SUM
      YCEN=YSUM/SUM
931   IF(IPRINT.EQ.0) GO TO 20
      WRITE(6,11) JH,JK,XMIN,NXM,NYM
11    FORMAT(//' FOR (',I4,',',I4,')  REFLECTION'/30X,'MINIMUM DENSITY
     1FROM TWO SIDES OF PERIMETER=',F6.1,' RASTER SIZE=',2I5)
      PERCNT=NABOVE*FACTOR
      WRITE(6,911)PERCNT,X
911   FORMAT(F8.1,' PERCENT OF RASTER ABOVE IMIN+',F6.2,
     1'*SQRT(IMIN)'/)
C     
      IXSTP=1
      IYSTP=1
      IF(NXM.LE.30)GO TO 1012
      IXSTP=2
      IYSTP=2
1012  DO 14 IY=1,NYM,IYSTP
14    WRITE(6,12)(JSTORE(IS,IX,IY),IX=1,NXM,IXSTP)
12    FORMAT(30I4)
      WRITE(6,31) XCEN,YCEN
31    FORMAT(70X,' REFINED POSN OF BACKGROUND CORR CENTRE IS',2F6.2)
C
      IF(ABNORM.LT.0.0)GO TO 1001
      JSTORE(IS,1,1)=-1
      IF (ISTOP.EQ.1)GO TO 120
      RETURN
C
1001  ABNORM=1.0
      JSTORE(IS,1,1)=-1
      IF(IPRINT.EQ.1)WRITE(6,1002)
1002  FORMAT(' REFLECTION EXPLICITLY EXCLUDED FROM LATTICE REFINE')
      IF(ISTOP.EQ.1)GO TO 120
      RETURN
C
950   ABNORM=1.0
      JSTORE(IS,1,1)=-1
      WRITE(6,951)JH,JK
951   FORMAT(' RASTER FOR REFLECTION',2I5,' GOES OUTSIDE MAX RAD')
      IF(ISTOP.EQ.1)GO TO 120
      RETURN
C
50    ABNORM=1.0
      PERCNT=(100.*NABOVE)/(NXM*NYM)
      JSTORE(IS,1,1)=-1
C   
      IF(IPRINT.EQ.1.OR.IPRINT.EQ.2) WRITE(6,51) PERCNT,X,JH,JK,
     1NXM,NYM
51    FORMAT(F8.1,' PERCENT OF RASTER ABOVE IMIN+',F6.2,
     1'*SQRT(IMIN)  CENTRE NOT CALC FOR REFLECTION',2I4,' RASTER',
     2' SIZE ',2I4)
      IF(IPRINT.NE.2)GO TO 1051
      IXSTP=1
      IYSTP=1
      IF(NXM.LE.30)GO TO 1013
      IXSTP=2
      IYSTP=2
1013  DO 1014 IY=1,NYM,IYSTP
1014  WRITE(6,12)(JSTORE(IS,IX,IY),IX=1,NXM,IXSTP)
1051  IF(ISTOP.EQ.1)GO TO 120
      RETURN
20    JSTORE(IS,1,1)=-1
      IF(ISTOP.EQ.1)GO TO 120
      RETURN
C
120   CONTINUE
C
      ODBTOT=0.0
      IOVERC=0
      DO 115 IX=1,NXM
      XCOORD=NXST+IX-1-X0
      XCOSQ=XCOORD**2
      DO 115 IY=1,NYM
      IYCURV=IY-1+NYST
      YCOORD=NYST+IY-1-Y0
      RAD=1.+SQRT(XCOSQ+YCOORD**2)
      IF(RAD.GT.(ROUT+1.))IOVERC=3
      IF(RAD.LE.RIN)IOVERC=3
      IRAD=RAD
      DRAD=RAD-IRAD
      ODB=ODBACK(IRAD)*(1-DRAD)+ODBACK(IRAD+1)*DRAD+YCURVE(IYCURV)
      JLOOK=ODB
      DJLOOK=ODB-JLOOK
CHENN>
C      ODBTOT=ODBTOT+LOOKUP(JLOOK)*(1-DJLOOK)+LOOKUP(JLOOK+1)*DJLOOK
C
C      write(6,'(''JLOOK1 = '',I8,2F15.3,I8)')JLOOK,ODB,DRAD,IRAD
C
      if(JLOOK.gt.0)then
        ODBTOT=ODBTOT+LOOKUP(JLOOK)*(1-DJLOOK)+LOOKUP(JLOOK+1)*DJLOOK
      else
        ODBTOT=ODBTOT+0.0
      endif
C
C      write(6,'(''JLOOK2 = '',I8,2F15.3,I8)')JLOOK,ODB,DRAD,IRAD
CHENN<
      IF(JLOOK.GT.IOVERL)IOVERC=1
      IF((JLOOK+1).LE.IUNDER)IOVERC=2
115   CONTINUE
      IF(IPRINT.EQ.0)GO TO 121
      WRITE(6,912)ODBTOT
912   FORMAT(60X,' SUMMATION OF RADIAL BACKGROUND CONTRIBUTION IS',
     1F10.0//)
C
121   JSTORE(IS,1,1)=-1
      RETURN
      END
C
C##############################################################################
C
      SUBROUTINE REFIND(X0,Y0,DX1,DY1,DX2,DY2,JH,JK,XREF,YREF,SINT,
     1NDATA,B3,XEWALD,YEWALD,X0OLD,Y0OLD,
     2DX1OLD,DY1OLD,DX2OLD,DY2OLD)
CHENN>
C      PARAMETER (NMAX=6000)
C
      PARAMETER (NMAX=48000)
CHENN<
      DIMENSION XREF(NMAX),YREF(NMAX),SINT(NMAX)
      INTEGER*2 JH(NMAX),JK(NMAX)
C
C     REFINE LATTICE PARAMETERS.
C
C
C
      DX1=DX1OLD
      DX2=DX2OLD
      DY1=DY1OLD
      DY2=DY2OLD
      X0=X0OLD
      Y0=Y0OLD
C
      DO 30 NAP=1,3
      S=0.
      SH=0.
      SK=0.
      SX=0.
      SY=0.
      SHH=0.
      SKK=0.
      SHK=0.
      SXH=0.
      SXK=0.
      SYH=0.
      SYK=0.
C
      DO 20 J=1,NDATA
      IF(XREF(J).EQ.0.0) GO TO 20
      RSQ=(JH(J)*DX1+JK(J)*DX2)**2 +
     1  (JH(J)*DY1+JK(J)*DY2)**2
      RADSQ=(XREF(J)-X0)**2 + (YREF(J)-Y0)**2
      XCORR=(XREF(J)-X0)*B3*RADSQ+XEWALD*RSQ
      YCORR=(YREF(J)-Y0)*B3*RADSQ+YEWALD*RSQ
      XOBS=XREF(J)-XCORR
      YOBS=YREF(J)-YCORR
      S=S+1.0
      SH=SH+JH(J)
      SK=SK+JK(J)
      SX=SX+XOBS
      SY=SY+YOBS
      SHH=SHH+JH(J)**2
      SKK=SKK+JK(J)**2
      SHK=SHK+JH(J)*JK(J)
      SXH=SXH+XOBS*JH(J)
      SXK=SXK+XOBS*JK(J)
      SYH=SYH+YOBS*JH(J)
      SYK=SYK+YOBS*JK(J)
20    CONTINUE
      IF (S.EQ.0)WRITE(6,920)
920   FORMAT(' NO REFLECTIONS TO USE IN REFIND')
      IF(S.EQ.0)STOP
      BOTTOM=(SH*SK-S*SHK)**2-(SK*SK-S*SKK)*(SH*SH-S*SHH)
      DX1=((SX*SK-S*SXK)*(SH*SK-S*SHK)-(SX*SH-S*SXH)*(SK*SK-S*SKK))/
     1 BOTTOM
      DX2=((SX*SH-S*SXH)*(SH*SK-S*SHK)-(SX*SK-S*SXK)*(SH*SH-S*SHH))/
     1 BOTTOM
      X0=(SX-SK*DX2-SH*DX1)/S
      DY1=((SY*SK-S*SYK)*(SH*SK-S*SHK)-(SY*SH-S*SYH)*(SK*SK-S*SKK))/
     1 BOTTOM
      DY2=((SY*SH-S*SYH)*(SH*SK-S*SHK)-(SY*SK-S*SYK)*(SH*SH-S*SHH))/
     1 BOTTOM
      Y0=(SY-SK*DY2-SH*DY1)/S
C
30    CONTINUE
      RETURN
      END
C
C##############################################################################
C
      SUBROUTINE SEARCH(X0,Y0,DX1,DY1,DX2,DY2,JH,JK,XREF,YREF,SINT,
     1NDATA,B3,XEWALD,YEWALD,X0OLD,Y0OLD,
     2DX1OLD,DY1OLD,DX2OLD,DY2OLD)
C
C     THIS SUBROUTINE WILL DO A SIMPLE SEARCH PROCEDURE ON B3 WITH
C     L.S. REFINEMENT OF X0,Y0,DX1,DY1,DX2,DY2 AND MINIMISATION OF L.
C     ALTERED SO THAT XD AND YD STAY FIXED
C     SEARCH B3 FROM 0 TO 2*B3INIT.
C     CHECK THAT THERE ARE ENOUGH REFLECTIONS TO USE.
C
CHENN>
C      PARAMETER (NMAX=6000)
C
      PARAMETER (NMAX=48000)
CHENN<
      INTEGER*2 JH(NMAX),JK(NMAX)
      DIMENSION XREF(NMAX),YREF(NMAX),SINT(NMAX)
      REAL LBEST,L
      INTEGER REDO
      LBEST=10000.
      B3INIT=0.00000001500
      B3STEP=0.00000000750
      NINT=0
      NCOORD=0
C      WRITE(6,911)XEWALD,YEWALD
C911   FORMAT(' XEWALD AND YEWALD',2F17.11)
      DO 10 J=1,NDATA
      IF(SINT(J).EQ.0.0) GO TO 5
      NINT=NINT+1
5     CONTINUE
      IF(XREF(J).EQ.0.0) GO TO 10
      NCOORD=NCOORD+1
10    CONTINUE
      WRITE(6,11) NINT,NCOORD
11    FORMAT(I6,' REFLECTIONS WITH MEASURED INTENSITY'/I6,
     1'  REFLECTIONS WITH ACCURATE CENTRES OF GRAVITY')
      IF(NCOORD.LE.10) GO TO 400
      WRITE(6,302)
      NCYCLE=0
19    NREDO=0 
20    REDO=0
      DO 100 IB3=1,5
      B3=B3INIT+(IB3-3)*B3STEP
      CALL REFIND(X0,Y0,DX1,DY1,DX2,DY2,JH,JK,XREF,YREF,SINT,
     1NDATA,B3,XEWALD,YEWALD,X0OLD,Y0OLD,
     2DX1OLD,DY1OLD,DX2OLD,DY2OLD)
C     WRITE(6,920)
C920   FORMAT(' AFTER REFIND')
C     WRITE(6,301)X0,Y0,DX1,DY1,DX2,DY2,XD,YD,B3
      L=0.
      N=0
      DO 50 I=1,NDATA
      IF(XREF(I).EQ.0.0) GO TO 50
      RADSQ=(XREF(I)-X0)**2+(YREF(I)-Y0)**2
      RSQ=(JH(I)*DX1+JK(I)*DX2)**2 + (JH(I)*DY1+JK(I)*DY2)**2
      XCALC=X0+JH(I)*DX1+JK(I)*DX2+B3*(XREF(I)-X0)*RADSQ+XEWALD*RSQ
      YCALC=Y0+JH(I)*DY1+JK(I)*DY2+B3*(YREF(I)-Y0)*RADSQ+YEWALD*RSQ
      A=(XREF(I)-XCALC)**2+(YREF(I)-YCALC)**2
      L=L+SQRT(A)
      N=N+1
50    CONTINUE
      L=L/N
      IF(L.GT.LBEST) GO TO 100
      LBEST=L
      NTEST=IB3
      B3BEST=B3
100   CONTINUE
      IF((NTEST.EQ.1).OR.(NTEST.EQ.5)) REDO=1
      CALL REFIND(X0,Y0,DX1,DY1,DX2,DY2,JH,JK,XREF,YREF,SINT,
     1NDATA,B3BEST,XEWALD,YEWALD,X0OLD,Y0OLD,
     2DX1OLD,DY1OLD,DX2OLD,DY2OLD)
      WRITE(6,301)X0,Y0,DX1,DY1,DX2,DY2,B3BEST,LBEST,
     1NCYCLE,NTEST
      B3INIT=B3BEST
      IF(REDO.NE.1) GO TO 51
      NREDO=NREDO+1
      IF(NREDO.GT.2)GO TO 51
      GO TO 20
51    NCYCLE=NCYCLE+1
      IF(NCYCLE.GT.5) GO TO 200
      B3STEP=B3STEP/4.
      GO TO 19
200   B3=B3BEST
      RETURN
400   WRITE(6,401)
401   FORMAT('  LESS THAN 10 REFLECTIONS FOUND-PARAMS MUST BE WRONG')
      STOP
301   FORMAT(2F10.2,4F10.3,F17.11,F8.4,I10,I3)
302   FORMAT('      X0        Y0       DX1       DY1       DX2',
     1'       DY2            B3        LAVER ',
     2'   NCYCLE IB3')
      END
C
C##############################################################################
C
	SUBROUTINE EMTILT(TL,TLTAXA,TLTANG,A0,B0,GAMMA0,A1,B1,GAMMA1)
C  CALCULATE TILT ANGLES FROM TILTED AND UNTILTED RECIPROCAL
C  CELL DIMENSIONS.
C
C  CONVENTION FOR MEASURING TILT AXIS TO ASTAR IS THAT THE ANGLE IS
C  FROM TILTAXIS TO ASTAR IN THE DIRECTION GIVEN BY ASTAR TO BSTAR
C  BEING POSITIVE.
C
	IMPLICIT REAL*8 (A-H,O-Z)
	REAL*4 TL,TLTAXA,TLTANG,A0,B0,GAMMA0,A1,B1,GAMMA1
	A=A0
	B=B0
	GAMMA=GAMMA0
	AT=A1
	BT=B1
	GAMMAT=GAMMA1
C
	COSG = DCOS(GAMMA*3.14159/180.0)
	COSGT= DCOS(GAMMAT*3.14159/180.0)
	SING = DSIN(GAMMA*3.14159/180.0)
C
	C1 = (A*A)/(AT*AT)
	C2 = (B*B)/(BT*BT)
	C3 = (A*B)/(AT*BT*COSGT)
	C4 = (A*B*COSG)/(AT*BT*COSGT)
C
	AX = C2*(C1/C3)*(C1/C3) - C1
	BX = C2 + 2.0*C2*(C1/C3)*((C1-C4)/C3) - C1
	CX = C2*((C1-C4)/C3)*((C1-C4)/C3)
	DISC = BX*BX - 4.0*AX*CX
	IF (DISC .LT. 0.0) GO TO 200
	PSQ1 = (-BX + DSQRT(DISC))/(2.0*AX)
	PSQ2 = (-BX - DSQRT(DISC))/(2.0*AX)
	IF (PSQ1 .LT. 0.0 .AND. PSQ2 .LT. 0.0) GO TO 210
	IF (PSQ1 .GT. 0.0) PP = DSQRT(PSQ1)
	IF (PSQ2 .GT. 0.0) PP = DSQRT(PSQ2)
	QQ = (C1/C3)*PP + ((C1-C4)/C3)*(1.0/PP)
	XK = DSQRT(C1*(PP*PP + 1.0))
C	WRITE(6,110) PP,QQ,XK
C110	FORMAT (' P,Q, AND SCALE FACTOR ',3F15.5)
	PHI = DATAN2(SING,QQ/PP-COSG)
	SINPHI = DSIN(PHI)
	COSPHI = DCOS(PHI)
	THETA = DATAN2(PP,SINPHI)
	TANTHE = DTAN(THETA)
	PHI = PHI*180.0/3.14159
	THETA = THETA*180.0/3.14159
C	WRITE(6,111) PHI,THETA
C111	FORMAT (' ANGLE FROM TILT AXIS TO ASTAR, AND TILT ANGLE',2F15.5)
C
C  CALCULATE ANGLE TO TILTED ASTAR ... THIS IS THE ANGLE ONE WOULD GET
C  DIRECTLY FROM THE FILM AND IS USED IN PICKTILT.
C
	COSPHI=COSPHI/(DSQRT(1.0+SINPHI*SINPHI*TANTHE*TANTHE))
	PHITLT=DACOS(COSPHI)*180.0/3.14159
C	WRITE(6,112) PHITLT
C112	FORMAT(' ANGLE FROM TILT AXIS TO TILTED ASTAR ON FILM',F15.5)
	TL=PHITLT
	TLTAXA=PHI
	TLTANG=THETA
	RETURN
200	WRITE (6,201)
201	FORMAT (' DISCRIMINANT LESS THAN ZERO')
	STOP
210	WRITE (6,211)
211	FORMAT(' TWO NEGATIVE ROOTS - SOMETHING WRONG')
	STOP
	END
C
C##############################################################################
C
      SUBROUTINE STATS(NDATA,JH,JK,ICORR,IOVER,
     .		X0,Y0,XCOORD,YCOORD,ASTAR,AHCORR,ROUT,RFA,NTOT)
      PARAMETER (MAXSLOTS=60)
      INTEGER*2 JH(1),JK(1),IOVER(1)
      DIMENSION ICORR(1)
      DIMENSION XCOORD(1),YCOORD(1)
      DIMENSION DSLOTC(MAXSLOTS),NSLOT(MAXSLOTS),TSLOTC(MAXSLOTS)
C
      DO 10 I=1,MAXSLOTS
      	NSLOT(I)=0
      	DSLOTC(I)=0.
      	TSLOTC(I)=0.
10    CONTINUE
      NREJECT=0
      DO 90 J=1,1+(2*NDATA/3)	! go well beyond centre
       DO 80 K=J+1,NDATA
        IF(JH(J).EQ.-JH(K)) THEN
      	  IF(JK(J).EQ.-JK(K)) THEN
      	    IF(IOVER(J).EQ.0.AND.IOVER(K).EQ.0) THEN
      		RADSQ=(XCOORD(J)-X0)**2+(YCOORD(J)-Y0)**2
      		ISLOT=0.999+(RADSQ/ROUT**2)*FLOAT(MAXSLOTS)
		IF(ISLOT.LT.1.OR.ISLOT.GT.6) STOP ' STATS resol wrong'
      		NSLOT(ISLOT)=NSLOT(ISLOT)+1
      		DSLOTC(ISLOT)=DSLOTC(ISLOT)+ABS(FLOAT(ICORR(J)-ICORR(K)))
      		TSLOTC(ISLOT)=TSLOTC(ISLOT)+
     .			AMAX1(1.,FLOAT(ICORR(J)+ICORR(K))/2.)
      	    ELSE
      		NREJECT=NREJECT+1
      	    ENDIF
      	  ENDIF
      	ENDIF
80     CONTINUE
90    CONTINUE
      WRITE(6,99)
99    FORMAT(//'      RESOLUTION DEPENDENT STATISTICS'/
     .	'      s   Resol    N   R_factor'/)
      NTOT=0
      DTOTC=0.
      TTOTC=0.
      DO 100 N=1,MAXSLOTS
       IF(NSLOT(N).GE.1) THEN
      	NTOT=NTOT+NSLOT(N)
      	DTOTC=DTOTC+DSLOTC(N)
      	TTOTC=TTOTC+TSLOTC(N)
      	R2=(AHCORR/(ASTAR*ROUT))*SQRT(FLOAT(MAXSLOTS)/FLOAT(N))
      	RFC=DSLOTC(N)/TSLOTC(N)
      	WRITE(6,101) N,R2,NSLOT(N),RFC
       ENDIF
100   CONTINUE
      IF(TTOTC.GT.0.0) RFA=DTOTC/TTOTC
      WRITE(6,102) NTOT,RFA,NREJECT
101   FORMAT(I7,F8.2,I5,F10.3)
102   FORMAT(/8X,'OVERALL',I5,F10.3,'  with',I5,
     .	' pairs over/underloaded '/)
      RETURN
      END
